Skip to content

Commit 17801d6

Browse files
odstrciltodstrcilt
authored
Update facit.py
1 parent 6538b52 commit 17801d6

File tree

1 file changed

+37
-12
lines changed

1 file changed

+37
-12
lines changed

‎aurora/facit.py‎

Lines changed: 37 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@ class FACIT:
4444
main ion density [:math:`1/m^3`]
4545
Nimp : nD array (...,nr) or float
4646
FSA impurity density [:math:`1/m^3`]
47-
Machi : nD array (...,nr) or float
47+
Mach_i : nD array (...,nr) or float
4848
main ion Mach number [-]
4949
Zeff: nD array (...,nr) or float
5050
effective charge [-]
@@ -62,10 +62,12 @@ class FACIT:
6262
major radius at magnetic axis [:math:`m`]
6363
qmag : 1D array (nr) or float
6464
safety factor [-]
65+
generalize to flux surface shaping as qmag_general = qmag * (r_min * B0 / d_torfluxov2pi_d_rmin)
66+
or equivalently qmag_general = r_min * B0 / d_polfluxov2pi_d_rmin
6567
rotation_model : int, optional
6668
if 0, use non-rotating limit from Fajardo PPCF (2022)
6769
if 1, use PS model from Maget PPCf (2020)
68-
if 2, use BP and PS model from Fajardo (subm. to PPCF)
70+
if 2, use BP and PS model from Fajardo PPCF (2023)
6971
Default is 0. 0 is recommended for low-Z impurities (e.g. B, N, Ne)
7072
Te_Ti : nD array (...,nr) or float, optional
7173
electron to ion temperature ratio [-]
@@ -184,7 +186,7 @@ def __init__(self, rho: (np.ndarray, float), \
184186
Zimp: (np.ndarray, float), Aimp: (int, float), \
185187
Zi: (np.ndarray, float), Ai: (int, float), \
186188
Ti: (np.ndarray, float), Ni: (np.ndarray, float), Nimp: (np.ndarray, float), \
187-
Machi: (np.ndarray, float), Zeff: (np.ndarray, float), \
189+
Mach_i: (np.ndarray, float), Zeff: (np.ndarray, float), \
188190
gradTi: (np.ndarray, float), gradNi: (np.ndarray, float), gradNimp: (np.ndarray, float), \
189191
invaspct: (float), B0: (float), R0: (float), qmag: (np.ndarray, float), \
190192

@@ -203,14 +205,16 @@ def __init__(self, rho: (np.ndarray, float), \
203205

204206
rho = np.maximum(1e-6, rho) # to avoid dividing by zero
205207

208+
assert np.size(invaspct) == 1 # new
209+
206210
if type(rho) is float or type(rho) is np.float64:
207211
nr = 1 # number of radial grid points
208212

209213
rho *= np.ones(nr)
210214
Ti *= np.ones(nr)
211215
Ni *= np.ones(nr)
212216
Nimp *= np.ones(nr)
213-
Machi *= np.ones(nr)
217+
Mach_i *= np.ones(nr)
214218
Zeff *= np.ones(nr)
215219
gradTi *= np.ones(nr)
216220
gradNi *= np.ones(nr)
@@ -252,9 +256,9 @@ def __init__(self, rho: (np.ndarray, float), \
252256
LnimpiNRL = 23. - np.log(Zi*Zimp*np.sqrt((Ni/1e6)*Zi**2+(Nimp/1e6)*Zimp**2)) + 1.5*np.log(Ti)
253257
LnimpimpNRL = 23. - np.log(Zimp**3*np.sqrt(2*(Nimp/1e6))) + 1.5*np.log(Ti)
254258

255-
# Braginskii collision times
259+
# Braginskii collision times (factor in front for NEO definition in contrast to NRL formulary definition)
256260

257-
Tauii = (self.eps_pi_fac*np.sqrt(mi)*Ti**1.5)/(Zi**4*Ni*LniiNRL) # ion-ion collision time [s]
261+
Tauii = ((1/(3.*np.sqrt(2.*np.pi)/4.))*self.eps_pi_fac*np.sqrt(mi)*Ti**1.5)/(Zi**4*Ni*LniiNRL) # ion-ion collision time [s]
258262
# set other collision times with respect to Tauii:
259263
Tauimpi = np.sqrt(Aimp/Ai)*((Zi**2*LniiNRL)/(Zimp**2*LnimpiNRL))*Tauii # impurity-ion collision time [s]
260264
Tauiimp = ((Zi**2*Ni*LniiNRL)/(Zimp**2*Nimp*LnimpiNRL))*Tauii # ion-impurity collision time [s]
@@ -273,7 +277,9 @@ def __init__(self, rho: (np.ndarray, float), \
273277

274278

275279
if rotation_model == 0:
276-
Machi *= 0.0
280+
Machi = 0.0
281+
else:
282+
Machi = Mach_i * 1.0
277283

278284
Machi2 = Machi**2 # auxiliary
279285

@@ -351,6 +357,7 @@ def __init__(self, rho: (np.ndarray, float), \
351357

352358

353359
e0imp = 1.0
360+
f_conv_fsa = 0.0
354361

355362

356363
elif rotation_model == 1:
@@ -409,7 +416,7 @@ def __init__(self, rho: (np.ndarray, float), \
409416
CclG = 1.0 + eps*deltan + 2*eps2
410417

411418
e0imp = 1.0
412-
419+
f_conv_fsa = 0.0
413420

414421

415422
elif rotation_model == 2:
@@ -421,8 +428,10 @@ def __init__(self, rho: (np.ndarray, float), \
421428

422429
if fsaout:
423430
e0imp = self.fluxavg(np.exp(Mzstar[...,None]**2*((RV**2 - (RV[:,0]**2)[:,None])/R0**2)), JV)
431+
f_conv_fsa = np.gradient(e0imp, amin*rho)/e0imp
424432
else:
425433
e0imp = np.ones(nr)
434+
f_conv_fsa = 0.0
426435

427436

428437

@@ -451,6 +460,11 @@ def __init__(self, rho: (np.ndarray, float), \
451460
self.Z = ZV
452461
self.J = JV
453462
self.B = BV
463+
self.e0imp = e0imp
464+
self.f_conv_fsa = f_conv_fsa
465+
self.g = g
466+
self.ft = ft
467+
self.eps = eps
454468

455469
# Pfirsch-Schlüter (PS)
456470

@@ -461,7 +475,7 @@ def __init__(self, rho: (np.ndarray, float), \
461475

462476
self.Vrz_PS = -self.Dz_PS*grad_ln_Nimp + self.Kz_PS*grad_ln_Ni + \
463477
self.Hz_PS*grad_ln_Ti # PS radial flux per particle [m/s]
464-
self.Vconv_PS = self.Kz_PS*grad_ln_Ni + self.Hz_PS*grad_ln_Ti # PS convective velocity [m/s]
478+
self.Vconv_PS = self.Kz_PS*grad_ln_Ni + self.Hz_PS*grad_ln_Ti + f_conv_fsa*self.Dz_PS # PS convective velocity [m/s]
465479

466480
# Banana-Plateau (BP)
467481

@@ -472,7 +486,7 @@ def __init__(self, rho: (np.ndarray, float), \
472486
self.Vrz_BP = -self.Dz_BP*grad_ln_Nimp + self.Kz_BP*grad_ln_Ni + \
473487
self.Hz_BP*grad_ln_Ti # BP radial flux per particle [m/s]
474488

475-
self.Vconv_BP = self.Kz_BP*grad_ln_Ni + self.Hz_BP*grad_ln_Ti # BP convective velocity [m/s]
489+
self.Vconv_BP = self.Kz_BP*grad_ln_Ni + self.Hz_BP*grad_ln_Ti + f_conv_fsa*self.Dz_BP # BP convective velocity [m/s]
476490

477491
# Classical
478492

@@ -482,7 +496,7 @@ def __init__(self, rho: (np.ndarray, float), \
482496

483497
self.Vrz_CL = -self.Dz_CL*grad_ln_Nimp + self.Kz_CL*grad_ln_Ni + \
484498
self.Hz_CL*grad_ln_Ti # CL radial flux per particle [m/s]
485-
self.Vconv_CL = self.Kz_CL*grad_ln_Ni + self.Hz_CL*grad_ln_Ti # CL convective velocity [m/s]
499+
self.Vconv_CL = self.Kz_CL*grad_ln_Ni + self.Hz_CL*grad_ln_Ti + f_conv_fsa*self.Dz_CL # CL convective velocity [m/s]
486500

487501

488502
# Total
@@ -491,11 +505,22 @@ def __init__(self, rho: (np.ndarray, float), \
491505
self.Kz = self.Kz_PS + self.Kz_BP + self.Kz_CL # coefficient of the main ion density gradient [m^2/s]
492506
self.Hz = self.Hz_PS + self.Hz_BP + self.Hz_CL # coefficient of the main ion temperature gradient [m^2/s]
493507

494-
self.Vconv = self.Kz*grad_ln_Ni + self.Hz*grad_ln_Ti # convective velocity [m/s]
508+
self.Vconv = self.Kz*grad_ln_Ni + self.Hz*grad_ln_Ti + f_conv_fsa*self.Dz# convective velocity [m/s]
495509
self.Vrz = self.Vrz_PS + self.Vrz_BP + self.Vrz_CL # radial flux per particle [m/s]
496510

497511
self.Flux_z = self.Vrz*Nimp # radial collisional flux [1/(m s)]
498512

513+
# Neoclassical
514+
515+
self.Dz_ncl = self.Dz_PS + self.Dz_BP # neoclassical diffusion coefficient [m^2/s]
516+
self.Kz_ncl = self.Kz_PS + self.Kz_BP # neoclassical of the main ion density gradient [m^2/s]
517+
self.Hz_ncl = self.Hz_PS + self.Hz_BP # neoclassical of the main ion temperature gradient [m^2/s]
518+
519+
self.Vconv_ncl = self.Kz_ncl*grad_ln_Ni + self.Hz_ncl*grad_ln_Ti + f_conv_fsa*self.Dz_ncl # convective velocity [m/s]
520+
self.Vrz_ncl = self.Vrz_PS + self.Vrz_BP # radial flux per particle [m/s]
521+
522+
self.Flux_z_ncl = self.Vrz_ncl*Nimp # radial neoclassical flux [1/(m s)]
523+
499524

500525
# asymmetries
501526

0 commit comments

Comments
 (0)