Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Strange behavior of TGLF with D-T ions #310

Open
pabloprf opened this issue Jul 25, 2023 · 37 comments
Open

Strange behavior of TGLF with D-T ions #310

pabloprf opened this issue Jul 25, 2023 · 37 comments

Comments

@pabloprf
Copy link
Contributor

When running TGLF for a plasma that has D and T as two separate ions, I realized that sometimes (like in the example I'll show below) the ion heat flux from one of the ions can be negative and almost of the same magnitude as the positive flux from the other ion. This produces weird results, particularly when scanning over one quantity (such as a/LTi) or when using in transport solvers (TGYRO/PORTALS).

The three input.tglf files I'll be discussing can be found here: https://www.dropbox.com/scl/fo/xym4q9hplxsr5ftjfj0ef/h?rlkey=fk4ybaneob1uzaxzffwr1653z&dl=0

run0 - D and T as separate species, produces the following output:

==========================================================
TGLF eb6b2025 [2023-07-13] OSX_MONTEREY Darwin x86_64
==========================================================
[Parsing data in input.tglf]
no mpi
  D(R) = -9.3057E-04  D(I) =  2.7701E-03
  kinetic species =          6  non-kinetic species =          0
      Gam/Gam_GB      Q/Q_GB  Q_low/Q_GB    Pi/Pi_GB      S/S_GB
elec  1.3399E-02  5.9458E-01  5.9458E-01  4.5486E-09  2.7920E-01
ion1 -4.8682E-03 -2.4302E+00 -2.4302E+00  1.5861E-03 -1.5602E-01
ion2  2.5543E-01  2.7434E+00  2.7434E+00 -2.4320E-03 -1.0958E-01
ion3 -1.5631E-02  3.8852E-03  3.8852E-03 -1.7702E-05  8.3690E-03
ion4 -7.4624E-05  1.2118E-04  1.2118E-04 -2.9375E-07  4.7091E-04
ion5 -4.6379E-02  1.1830E-01  1.1830E-01 -6.5009E-05 -1.9517E-02

run1 - Changing the masses so that both ions have a mass of A=2.5 (1.25 in normalized units), produces:

[Parsing data in input.tglf]
no mpi
  D(R) = -9.3057E-04  D(I) =  2.7701E-03
  kinetic species =          6  non-kinetic species =          0
      Gam/Gam_GB      Q/Q_GB  Q_low/Q_GB    Pi/Pi_GB      S/S_GB
elec  1.4455E-01  5.9232E-01  5.9232E-01 -1.1315E-08  2.7758E-01
ion1  9.7480E-02  5.8076E-01  5.8076E-01  2.4496E-05 -1.2957E-01
ion2  9.7480E-02  5.8076E-01  5.8076E-01  2.4496E-05 -1.2957E-01
ion3 -5.2902E-03 -1.9522E-03 -1.9522E-03  1.3562E-05  6.3124E-03
ion4 -3.0572E-05  1.5666E-05  1.5666E-05  9.9062E-09  2.6232E-04
ion5 -6.3447E-04  5.9251E-02  5.9251E-02  3.1114E-05 -2.1984E-02

When we compare run0 (blue) and run1 (red), we can see that the issue is a few low-k modes that become stable in run0 and that lead to a negative portion of the ion heat flux spectrum.
image

This was run with SAT2 but I don't think that should matter, since the differences are observed in the linear stability plots already.

This can have a quite strong effect on parameter scans, as seen here for a a/LTi (changing all ion species simultaneously):
image

Does someone has an idea of why this is happening? As far as I can tell, the parameters of this run are not out of the ordinary.
This may be of interest to @cholland @mmuraca @AudreySaltzman @nthoward
@pabloprf
Copy link
Contributor Author

In case someone is wondering, having a single species with mass A=2.5 and double the concentration (instead of two species with A=2.5, as in run1) produces somewhat similar results to run1:

run2

[Parsing data in input.tglf]
no mpi
  D(R) = -9.3057E-04  D(I) =  2.7701E-03
  kinetic species =          5  non-kinetic species =          0
      Gam/Gam_GB      Q/Q_GB  Q_low/Q_GB    Pi/Pi_GB      S/S_GB
elec  1.2525E-01  6.3724E-01  6.3724E-01  1.0573E-09  2.9803E-01
ion1  2.3485E-01  1.0200E+00  1.0200E+00  3.0694E-05 -2.7995E-01
ion2 -9.3172E-03 -5.2380E-03 -5.2380E-03  1.4725E-05  7.6465E-03
ion3 -5.5534E-05  1.7702E-05  1.7702E-05  7.6100E-08  3.1176E-04
ion4 -1.1485E-02  5.2460E-02  5.2460E-02  3.2914E-05 -2.2920E-02

@jcandy
Copy link
Member

jcandy commented Jul 26, 2023

Interesting. @gmstaebler needs to weigh in here.

Something very strange is happening to TGLF when ions have similar but unequal mass. Consider this input file:

# Control parameters
# ------------------

USE_TRANSPORT_MODEL     = True
GEOMETRY_FLAG           = 1
KY                      = 0.3
SAT_RULE                = 2
USE_BPER                = True
UNITS                   = CGYRO

# Geometry parameters
# ------------------

RMIN_LOC                = 0.55
RMAJ_LOC                = 3.3
DRMAJDX_LOC             = -0.07
Q_LOC                   = 1.03
KAPPA_LOC               = 1.4
P_PRIME_LOC             = -0.0006
Q_PRIME_LOC             = 2.4

# Plasma parameters
# ------------------

NS                      = 3
BETAE                   = 0.003
XNUE                    = 0.02

# Species
# -------

# Specie #1 (electrons)
ZS_1         = -1.0
MASS_1       = 0.000272445
RLNS_1       = 0.12
RLTS_1       = 0.6
TAUS_1       = 1.0
AS_1         = 1.0

# Specie #2 (thermal ion)
ZS_2         = 1.0
MASS_2       = 1.01
RLNS_2       = 0.12
RLTS_2       = 1.35
TAUS_2       = 1.0
AS_2         = 0.5

# Specie #3 (thermal ion)
ZS_3         = 1.0
MASS_3       = 1.0
RLNS_3       = 0.12
RLTS_3       = 1.35
TAUS_3       = 1.0
AS_3         = 0.5

Note that the ion masses are 1.0 and 1.01. This gives completely absurd fluxes:

      Gam/Gam_GB      Q/Q_GB  Q_low/Q_GB    Pi/Pi_GB      S/S_GB
elec  3.4859E-02  6.1907E-01  6.1907E-01 -1.8688E-07  3.6909E-01
ion1  7.5021E+00 -1.2530E+02 -1.2530E+02  8.1175E-03 -7.9478E-01
ion2 -7.4673E+00  1.2675E+02  1.2675E+02 -7.8917E-03  4.3071E-01

Making only a tiny change to equal mass: 1.0 and 1.0 gives a sensible result:

    Gam/Gam_GB      Q/Q_GB  Q_low/Q_GB    Pi/Pi_GB      S/S_GB
elec  1.0600E-01  5.9479E-01  5.9479E-01 -3.3726E-08  3.6600E-01
ion1  5.2999E-02  6.6853E-01  6.6853E-01  6.3427E-05 -1.8051E-01
ion2  5.2999E-02  6.6853E-01  6.6853E-01  6.3427E-05 -1.8051E-01

@gmstaebler
Copy link
Member

gmstaebler commented Jul 26, 2023 via email •
edited by jcandy
Loading

@pabloprf pabloprf reopened this Jul 26, 2023
@pabloprf
Copy link
Contributor Author

(sorry I clicked "send and close" before I could finish writing)

Hi @gmstaebler, thanks for the discussion!

I'm looking at the simple case that @jcandy proposed (I also run it electrostatic, to make it simpler) and I found the following. The linear growth rate and real frequency spectra are almost the same for both cases (blue with the unequal masses, red equal masses):
image
Even though the growth rates and real frequency spectra are basically the same, the wavefunction is not. I'm plotting here at ky~0.67 (second to highest ky mode):
image
And this different wavefunction structure has a strong effect on the Qi spectra. In blue, the two ions have very large and opposite in sign values:
image
While zooming-in the red case has a normal-looking spectrum:
image

If it's a compilation issue, we are facing the same problem at MIT and GA.

@gmstaebler
Copy link
Member

gmstaebler commented Jul 26, 2023 via email

@pabloprf
Copy link
Contributor Author

This is the result of a scan of MASS_2 from 1.0 to 1.5. In the following plots, the color transitions from red (1.0) to blue (1.5).
image

There are significant discontinuities in the total flux, which happens because a mode at ky~0.2 gets stabilized when the mass of the ion is (interestingly) above 1.25. The contributions to the changes in the total flux come from changes in the Qi spectrum below ky<0.4.

However, this discontinuity because of the ky=0.2 is a separate issue. At higher ky (0.6-0.8) is where we observe the strange behavior of the large-magnitude, opposite-sign contributions from both ions to the Qi spectrum. In the right most plot I have taken the minimum of the first ion flux spectrum and the maximum of the second ion flux spectrum, to evaluate if this behavior becomes stronger as the masses are more equal to each other. Indeed that's the case, as the mass approaches 1.0, the fluxes become asymptotically(?) very large.

@gmstaebler
Copy link
Member

gmstaebler commented Jul 27, 2023 via email

@pabloprf
Copy link
Contributor Author

HI @gmstaebler,

@nthoward and I have been running hundreds of nonlinear CGYRO simulations with ions of similar mass (D-T) and I don't think we have never seen a negative ion heat flux spectrum for any of the ions...

In any case, another piece of relevant information to explore why this is happening is that it's related to the saturation rule and not to the linear physics. I run the exact same case with SAT0 (doing USE_PRESETS = .FALSE. in tglf_startup.f90 so that the variables XNU_MODEL, WDIA_TRAPPED and UNITS are the same between SAT0 and SAT2) and the results indicate that the linear growth rates, frequencies and wavefunctions are exactly the same between the two (as expected), yet the odd negative-positive Qi spectrum for ion1 and ion2 only appears in SAT2.

[Blue is SAT2 and red is SAT0]
Growth rates and real frequencies overlap:
image
Wavefuntions overlap:
image
Qi spectrum for SAT2 has the weird behavior:
image
Zoom shows SAT0 does not have that behavior:
image

@gmstaebler
Copy link
Member

gmstaebler commented Jul 31, 2023 via email

@pabloprf
Copy link
Contributor Author

I'm using a similar input file that Jeff sent previously:

#-------------------------------------------------------------------------
# TGLF input file
#-------------------------------------------------------------------------

# Control parameters
# ------------------

USE_TRANSPORT_MODEL     = True
GEOMETRY_FLAG           = 1
KY                      = 0.3
SAT_RULE                = 2
USE_BPER                = False
UNITS                   = CGYRO
XNU_MODEL               = 3
WDIA_TRAPPED            = 1.0


# Geometry parameters
# ------------------

RMIN_LOC                = 0.55
RMAJ_LOC                = 3.3
DRMAJDX_LOC             = -0.07
Q_LOC                   = 1.03
KAPPA_LOC               = 1.4
P_PRIME_LOC             = -0.0006
Q_PRIME_LOC             = 2.4


# Plasma parameters
# ------------------

NS                      = 3
BETAE                   = 0.003
XNUE                    = 0.02


# Species
# -------

# Specie #1 (electrons)
ZS_1         = -1.0
MASS_1       = 0.000272445
RLNS_1       = 0.12
RLTS_1       = 0.6
TAUS_1       = 1.0
AS_1         = 1.0
VPAR_1       = 0.0
VPAR_SHEAR_1 = 0.0
VNS_SHEAR_1  = 0.0
VTS_SHEAR_1  = 0.0

# Specie #2 (thermal ion)
ZS_2         = 1.0
MASS_2       = 1.01
RLNS_2       = 0.12
RLTS_2       = 1.35
TAUS_2       = 1.0
AS_2         = 0.5
VPAR_2       = 0.0
VPAR_SHEAR_2 = 0.0
VNS_SHEAR_2  = 0.0
VTS_SHEAR_2  = 0.0

# Specie #3 (thermal ion)
ZS_3         = 1.0
MASS_3       = 1.0
RLNS_3       = 0.12
RLTS_3       = 1.35
TAUS_3       = 1.0
AS_3         = 0.5
VPAR_3       = 0.0
VPAR_SHEAR_3 = 0.0
VNS_SHEAR_3  = 0.0
VTS_SHEAR_3  = 0.0

@pabloprf
Copy link
Contributor Author

pabloprf commented Aug 3, 2023

Some more info on this issue.

The QL flux spectrum (top row) is exactly the same for SAT0 and SAT2. Both show the very odd negative-positive Qi for the two ion heat fluxes. As expected, as the linear physics was the same.
However, the reason why that doesn't show up in the Sum Flux spectrum (bottom row) for SAT0 is because the potential spectrum turns out to be equal to zero at that specific ky mode in SAT0. It's not zero in SAT2, which ends up amplifying the QL spectrum observation.

image

This, to me, brings up back to the linear solver as the source for this behavior. It was just that SAT0 made the issue disappear because the saturated potential at that wavenumber was exactly zero.

Does this make sense? Am I interpreting correctly the content of the QL_flux_spectrum and sum_flux_spectrum output files?

@mmuraca
Copy link

mmuraca commented Aug 3, 2023 via email

@pabloprf
Copy link
Contributor Author

pabloprf commented Aug 3, 2023

Marco, delta phi for SAT0 is zero at the wavenumber with the very large QL Qi flux: ky~0.64

@gmstaebler
Copy link
Member

gmstaebler commented Aug 3, 2023 via email

@pabloprf
Copy link
Contributor Author

pabloprf commented Aug 3, 2023

Hi Gary,
Actually, the growth rate is not zero at ky~0.64. It's only that the potential intensity in SAT0 is zero. In SAT2 isn't
This is a more complete plot with the growth rate spectrum on the top left and moving the potential spectrum in the bottom left. Vertical line is the mode that has the large QL Qi fluxes of opposite sign for each ion:
image

@tomneiser
Copy link
Contributor

Hi Pablo, these are interesting findings - thanks for sharing. I was curious about Gary's question regarding linear CGYRO, and also if this is exclusive to TGLF SAT2 default settings. I have attached an example that shows that TGLF SAT1 default settings can have the same behavior as the TGLF SAT2 defaults (see low-ky vertical line in my plot). CGYRO also has the same behavior at high-ky, but after more testing I found this behavior generally only occurs when the linear mode is stable and so does not affect QLGYRO fluxes (corroborating Pablo and Nathan's nonlinear results).

DT_TGLF_SAT1_vs_linCGYRO

@pabloprf
Copy link
Contributor Author

Interesting, thanks @tomneiser!
So if I understood correctly your plot, you verified that the strange negative-positive low-k Qi QL fluxes that happen in TGLF do not occur with linear CGYRO. And that they are independent from the setup of the TGLF linear solver (SAT1 vs SAT2 defaults).
Do you have an idea of what we should investigate next to resolve this? Could this be physical (linear GF equations exhibiting this vs linear GK equations do not) or a bug somewhere?

@jcandy
Copy link
Member

jcandy commented Aug 14, 2023

When I first saw the onset of anomalous behaviour that got worse as m1-m2 -> 0, I thought it might be related to the TGLF eigensolver having problems when two eigenvalues are nearly equal

@gmstaebler
Copy link
Member

gmstaebler commented Aug 14, 2023 via email •
edited by jcandy
Loading

@jcandy
Copy link
Member

jcandy commented Aug 14, 2023

@gmstaebler this is something different, since are there are not two nearby physical eigenvalues in this case. The eigenvalue and eigenmodes should be smooth functions of m2-m1. However, TGLF seems to have an unphysical sensitivity to m2-m1. Of course it works properly in the limit m2-m1=0, so that is a clue.

Also, you can reply to these issues directly using github (the interface is nice):

#310

@tomneiser
Copy link
Contributor

@pabloprf The blue and green lines in the left panel of my plot are quasi-linear weights calculated with linear CGYRO, and they do show this same negative-positive QL weights behavior as TGLF. However, whenever I see this behavior for CGYRO the linear mode is stable (I set the linear CGYRO growth rate to zero whenever gamma<0, e.g. blue line in the right panel). And yes, the strange behavior happens independent of the TGLF linear solver.

I hope it's a helpful clue that this behavior also happens in linear CGYRO (but is harmless due to zero growth rates). The linear TGLF growth rate spectrum of the most unstable mode (green line in the right panel) shows a branch jump at the offending ky, so it could be a weakness of the eigenvalue solver when there are two nearby modes. @pabloprf's most recent detailed plot (top left panel) shows two modes with similar eigenvalues in the offending region, so there could be some mass-dependent coupling between nearby modes in TGLF that's causing this. To find the bug I would start by looking at mass-dependent terms in the QL weight calculation

@pabloprf
Copy link
Contributor Author

An update on my investigation of this this.
The issue that happens at ky~0.645 appears already in the eigenvector solution from the linear solver, and it's made more evident when comparing the energy_weight from

gacode/tglf/src/tglf_LS.f90

Lines 803 to 804 in 08ff9dc

energy_weight(is,1) = energy_weight(is,1) &
+ REAL(xi*CONJG(phi(i))*p_tot(is,i))

and the phi_norm from

gacode/tglf/src/tglf_LS.f90

Line 712 in 08ff9dc

phi_norm = phi_norm + REAL(phi(i)*CONJG(phi(i)))

The magnitude of energy_weight for the ions is much larger than phi_norm, which explodes in the construction of the QL weight here:

gacode/tglf/src/tglf_LS.f90

Line 843 in 08ff9dc

energy_weight(is,j) = as(is)*taus(is)*1.5*ky*energy_weight(is,j)/phi_norm

In fact, this behavior is already evident when plotting the most unstable eigenmode from the case with equal masses (blue) and the case with 1% difference (red). Note the 6 orders of magnitude difference in scale between the two plots, even if both correspond to the same eigenvalue (gamma~0.0096):

image

image

If I set a minimum value of phi_norm to 1.E-9 instead of 1.E-12 here:

gacode/tglf/src/tglf_LS.f90

Line 619 in 08ff9dc

epsilon1 = 1.E-12

Then the issue at that wavenumber disappears and the case with MASS_3=1.01 gets very similar to MASS_3=1.0

I hope this is relevant for the experts in the room to figure out a solution!

@pabloprf
Copy link
Contributor Author

(To give some numbers to the magnitudes of energy_weight and phi_norm that I mentioned on my previous comment)

I added these before line 841 in tglf_LS.f90:

if(is.eq.1)write(*,*)"ky = ",ky,"energy_weight(1)",energy_weight(1,1),"phi_norm",phi_norm
if(is.eq.2)write(*,*)"ky = ",ky,"energy_weight(2)",energy_weight(2,1),"phi_norm",phi_norm

And for the most unstable mode it results:

 ky =   0.64516129032258074      energy_weight(1)   3.7178003918532221E-013 phi_norm   3.7260760284206491E-012
 ky =   0.64516129032258074      energy_weight(2)   7.0467907248049583E-009 phi_norm   3.7260760284206491E-012

The large magnitude of energy_weight(2) compared to the small magnitude of phi_norm is what's causing the trouble at that wavenumber, and the reason why it gets resolved by giving a minimum value to phi_norm through epsilon1.

@jcandy
Copy link
Member

jcandy commented Aug 15, 2023

Ah, so epsilon1 is what corrects the result in the limit m2-m1->0. I was puzzled at why the problem suddenly disappeared as the difference decreased.

@pabloprf
Copy link
Contributor Author

Without the epsilon1 limit, the problem because more pronounced as m2-m1->0 because phi_norm decreases more strongly as than energy_weight(2) , but actually at some point it must be reversed because for m2=m1, even without the epsilon1 limit, I think the problem is solved by itself:

MASS_3=1.01

 B, ky =   0.64516129032258074      energy_weight(1)   3.7178003918532221E-013 phi_norm   3.7260760284206491E-012
 B, ky =   0.64516129032258074      energy_weight(2)   7.0467907248049583E-009 phi_norm   3.7260760284206491E-012

MASS_3=1.005

 B, ky =   0.64516129032258074      energy_weight(1)   2.3577462291105267E-014 phi_norm   2.3598653379161407E-013
 B, ky =   0.64516129032258074      energy_weight(2)   8.8498295403213366E-010 phi_norm   2.3598653379161407E-013

MASS_3=1.0025

 B, ky =   0.64516129032258074      energy_weight(1)   1.4944045232664915E-015 phi_norm   1.4947886729277071E-014
 B, ky =   0.64516129032258074      energy_weight(2)   1.1163277323531631E-010 phi_norm   1.4947886729277071E-014

MASS_3=1.001

 B, ky =   0.64516129032258074      energy_weight(1)   3.8961813810949645E-017 phi_norm   3.8963366655651828E-016
 B, ky =   0.64516129032258074      energy_weight(2)   7.2558029665949580E-012 phi_norm   3.8963366655651828E-016

MASS_3=1.0

 B, ky =   0.64516129032258074      energy_weight(1)  -7.0899454289492581E-025 phi_norm   1.5149580944418587E-023
 B, ky =   0.64516129032258074      energy_weight(2)   6.4645565536917493E-024 phi_norm   1.5149580944418587E-023

@pabloprf
Copy link
Contributor Author

This looks to me a numerical problem of the eigenvalue solver that could be mitigated by increasing epsilon1 (do people agree?). However, how much to increase it, I don't know. Maybe @gmstaebler can comment on how epsilon1=1.E-12 was chosen or what determines it.
(btw, as Jeff mentioned, the github interface to see these messages is significantly better than email)

@gmstaebler
Copy link
Member

gmstaebler commented Aug 16, 2023 via email •
edited by jcandy
Loading

@gmstaebler
Copy link
Member

gmstaebler commented Aug 16, 2023 via email •
edited by jcandy
Loading

@tomneiser
Copy link
Contributor

This is helpful information and good news, thank you Gary!

@pabloprf
Copy link
Contributor Author

This is good! Thanks Gary.
I don't know what the GFS solver is though. Is that available to run TGLF with? What changed?

@gmstaebler
Copy link
Member

gmstaebler commented Aug 16, 2023 via email •
edited by tomneiser
Loading

@gmstaebler
Copy link
Member

Pablo,Here is what I get running the stand alone TGLF for this case
TGLF 94e45ce [2023-07-05] OSX_VENTURA Darwin arm64

[Parsing data in input.tglf]
no mpi
D(R) = 2.5817E-03 D(I) = 4.7561E-03
kinetic species = 3 non-kinetic species = 0
Gam/Gam_GB Q/Q_GB Q_low/Q_GB Pi/Pi_GB S/S_GB
elec 4.6560E-02 6.4743E-01 6.4743E-01 -1.6736E-07 3.9963E-01
ion1 5.1592E+00 -1.0381E+02 -1.0381E+02 9.2174E-03 -7.2527E-01
ion2 -5.1126E+00 1.0561E+02 1.0561E+02 -9.0987E-03 3.2564E-01

@gmstaebler
Copy link
Member

Note that SAT0 uses the total eigenvector as a normalization for the QL weights not the potential so that could make a big difference. The SAT2 linear eigenmodes are also different but SAT1 and SAT0 should be the same

@pabloprf
Copy link
Contributor Author

Hi Gary,
Thanks for the info about what GFS is! looking forward to seeing the publication. Is the plan to have it implemented as the default eigensolver in TGLF? If so, do you have a plan on timeline for that?

Those values of running TGLF standalone are a bit different from what Jeff sent earlier in this thread, but your values also show the positive-negative Qi pair, which give very different results when compared to the GFS results you showed earlier.

@pabloprf
Copy link
Contributor Author

Also, while GFS is not available as a solver in TGLF, what do you recommend to avoid the potential bug?

@gmstaebler
Copy link
Member

gmstaebler commented Aug 16, 2023 via email

@gmstaebler
Copy link
Member

gmstaebler commented Aug 17, 2023 via email •
edited by tomneiser
Loading

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
5 participants
@jcandy @gmstaebler @tomneiser @pabloprf @mmuraca