The NCLASS_PT stand-alone driver supplied with the NCLASS module provides sample input for a single call to NCLASS (a single flux surface). It uses a simplified large aspect ratio elongated plasma approximation to generate the required geometric input. It then gives examples of reconstructing the fluxes, bootstrap current, and rotation velocities from the NCLASS module output and writes a summary of the results of the test case. If the user is interested in incorporating NCLASS as a module in another code, the driverroutine provides valuable examples of how to manipulate the output to generate other related quantities.
For inclusion of the NCLASS module in a transport or other analysis code, there are additional modules available from the author for reading VMEC or EFIT output files, and others for calculating the required flux surface averaged geometric quantities from either of those representations.
Significant effort has gone into making the NCLASS module computationally efficient so it can be used for routine analysis in transport codes. It has been incorporated into several transport codes (e.g., TRANSP, ONETWO, ASTRA) and also used extensively in another stand-alone form (FORCEBAL code) for analysis of data (e.g., at PPPL, GA, JET, TEXTOR, Tore Supra). Although the numerical procedures are rather sophisticated, it is computationally robust.
All units are in mks with temperatures in keV. The radial variable can be in any unit of the user's choice as long as all input to NCLASS is consistent (see Radial grid). There is no way to check this consistency inside NCLASS. Many of the output quantities also carry the arbitrary unit of the radial variable (see NCLASS Variables or the comment blocks at the beginning of each routine).
File names and routines are rendered as bold uppercase in a
typewriter font: FILENAME.EXT
Coding and variable names are rendered as bold lowercase in a typewriter
font: p_ft
Greek characters appear as (if the Symbol font is available these are greek letters): abcd
README_NCLASS_PT.HTML
:
README_NCLASS_PT.HTML | |--NCLASS_PT_TOP.HTML | |--NCLASS_PT_INFO.HTML | |--NCLASS_PT_INDEX.HTML | |--COMPUTER.GIF | |--COUNTER.GIFThe GIF files have been added for entertainment, courtesy of my daughter, Michele, who convinced me HTML was a snap to learn. In fact, it is almost as easy as LaTeX.
NCLASS_MOD.FORthe driver,
NCLASS_PT_DR.FORthree PARAMETER 'include' files,
PAMX_MI.INC PAMX_MS.INC PAMX_MZ.INCa few auxiliary function routines,
U_ERF.FOR U_LU_BACKSUB.FOR U_LU_DECOMP.FORarray manipulation routines,
RARRAY_COPY.FOR RARRAY_SUM.FOR RARRAY_ZERO.FORand a module that control the 'write' formats,
WRITE_MOD.FOREach routine has a description block at the beginning that covers functionality, and all I/O. A flow and dependency diagram for NCLASS_PT is shown under Flow Diagram
The parameters have been set up to handle a maximum of nine isotopes at a time, a maximum of forty species (populated charge states of all isotopes) and Z up to 18. Although this is more than adequate for most problems, NCLASS has been applied to problems that exceed each of these limits. A diagnostic error is generated by NCLASS if the problem calls for exceeding any of these limits.
MAKEFILE
IN_NCLASS_PT.NMLAll input variables are set to the test case values as defaults in the driver. Below is a description of the NAMELIST variables, and whether they are used directly for input to NCLASS (see NCLASS input) or used by the driver routine (see Driver input) to generate other required NCLASS input, control the output, etc.
k_out-option for output to nout (-) =1 errors only =2 errors and results =else no output p_eps-inverse aspect ratio (-) p_q-safety factor (-) r0-major radius (m) a0-minor radius, scale length for rho (m) e0-axial elongation (-) bt0-axial toroidal field (T) q0-axial safety factor (-)
k_order-order of v moments to be solved (-) =2 u and q =3 u, q, and u2 =else error k_potato-option to include potato orbits (-) =0 off =else on c_den-density cutoff below which species is ignored (/m3) p_grphi-radial electric field Phi' (V/rho) p_gr2phi-radial electric field gradient Psi'(Phi'/Psi')' (V/rho2) amu_i(i)-atomic mass number of i (-) temp_i(i)-temperature of i (keV) den_iz(i,z)-density of i,z (/m3)The
p_*
values are for the flux surface of
interest, i.e., radial variables when NCLASS is called in a loop
over radial nodes. (For a discussion of naming conventions, see
Name Conventions). The *0
values are global or axial values to generate the potato orbit
parameters and other approximate geometric quantities.
a0
is used to generate characteristic gradient
scales. Phi
denotes the electric potential. The radial
electric field and its derivative are expressed in terms of the
electric potential because NCLASS uses flux surface averaged
quantities in its equations. Psi
is the poloidal
flux contained within the flux surface of interest (not divided
by 2p).
NCLASS_PT sets up the gradients and metrics using a radial
grid, r, that is proportional to the square root
of the toroidal flux and normalized to a0
so that it
has units of meters, and '
designates derivatives with
respect to r.
amu_i
and temp_i
input, and the
isotopic populations by den_iz
.
amu_i
< 0.5) to calculate the
transport resulting from ion-impurity collisions only. The
resulting fluxes should be ambipolar. This case can be run by
modifying the input in the IN_NCLASS_PT.NML
file.
amu_i(2)
= amu_i(3)
= 2.0 with
different temperatures and densities for isotopes 2 and 3
governed by den_iz
and temp_i
input).
This case can be run by modifying the input in the
IN_NCLASS_PT.NML
file.
k_potato
= 1. This provides an empirical
modification of the viscosity near the magnetic axis to match
theoretical limits. With large potato orbit corrections, the
trapped fraction does not vanish at the axis. Although the
potato orbit results have appeared reasonable in various tests
and applications to TFTR, DIII-D and ITER the viscosity has not
been benchmarked against more detailed kinetic calculations to
verify either the radial extent of the orbit corrections or the
form.
p_fm
to zero.
p_ft
) to
zero.
p_fm
and p_ft
to zero.
p_fhat
. (dY/dr
enters through the denominator of the NCLASS input variable
p_fhat
). For other grid representations (e.g.,
where r is proportional to the square root of
the toroidal flux as in NCLASS_PT), avoid calling NCLASS at
the origin where the variable p_fhat
is singular.
The simple geometry and choice of grid in NCLASS_PT do not
allow the user to calculate axial modifications without
revising the driver.
OUT_NCLASS_PT.DATIf
k_out
= 1 is selected, the only output will be
warnings and error messages. For k_out
= 2,
warnings and error messages plus results will be output. Other
choices of k_out
eliminate all output.
k_potato
.
The parallel current response to the external force has been
added. An example of the use of this feature is that it
gives the response current to neutral beam current drive when
the beam parallel forces are input through fex_iz
.
The calculations of one of the internal routines that had minimal functionality have been absorbed by the interface routine.
All I/O for the NCLASS module is now passed as arguments through the interface routine. This facilitated the removal of all COMMON blocks.
k_order
= 2 is specified,
reducing the computational time by nearly a factor of two. For
most applications only two moments are necessary. In
collisional plasmas (Pfirsch-Schlüter regime) it can be
important to include the third moment terms. Not all the
important third moment terms have been implemented in NCLASS
(see closure discussion in the Hirshman and Sigmar review, it's
messy). The retention of the third order terms is to accomodate
further development, primarily for edge plasmas.
Some of the internal variables have been redefined as 'species' arrays to reduce the checking for significantly populated states. The mapping from (mass, charge state) to populated species is done in the interface, used extensively internally, and then passed out to help the user interpret the results. Because the (mass,charge state) arrays are typically very sparse (usually only a single charge state is used for each isotope) this eliminates a lot of unnecessary storage, time consumption in passing arrays through arguments, and zeroing out large arrays.
All variables are declared and identified as either input, output, or local in all routines in both NCLASS and NCLASS_PT.
All COMMON blocks have been removed.
Error checking of input has been provided so that array limits governed by PARAMETER statements are not exceeded.
Error messages are all passed by a flag from NCLASS to the driver routine.
All WRITE formats are isolated in separate routines and
controlled by a user option
k_out
to NCLASS_PT.
Extensive use has been made of variable and routine naming conventions (see Name Conventions).
Although it is not possibe to adhere to the goal of dynamic allocation of variables in FORTRAN 77, NCLASS and NCLASS_PT do make use of PARAMETER statements to allocate the storage needed for all major arrays. In NCLASS this is a fairly important consideration because of the potentially large arrays when multiple charge states of several impurities are of interest. These limits are checked during each calculation and error messages are generated if the calculation will exceed a limiting parameter.
capn_ii(3,3,i1,i2)
, was corrected to obey the required
symmetry constraint for unequal temperatures. The expression is
correct in the primary NCLASS reference (see
NCLASS equations).
iflag
is set to a negative
value to warn the user in case the null or unity input values were
input by mistake (see NCLASS output).
Such limits include calculations at the axis (no trapped particles,
p_ft=0
) or at the edge of a unity aspect ratio
plasma (no circulating particles, p_ft=1.0
). The
banana contributions can be turned off at any radius with
p_ft=0
to facilitate a calculation with only
Pfirsch-Schlüter contributions to the viscosity. The
Pfirsch-Schlüter contribution to viscosity can be turned off
by setting all three poloidal moments of p_fm(j)=0.0
.
If both p_ft=0
and p_fm(j)=0.0
, the
Spitzer problem is solved.
The viscosities are calculated by numerically integrating over velocity space and are continuous over all collisionality regimes and aspect ratios. In that regard the model is more advanced than the one presented in Hirshman and Sigmar's review paper and other codes that use the matrix approach for calculating neoclassical properties (e.g., work by Tani, et al., Kessel, Ernst and a few others). Other advanced features include banana orbit squeezing, potato orbit effects and additional force terms to accomodate neutral beam, charge exchange, toroidal field ripple and anomalous toroidal drag forces. Nonetheless, the model must still be considered only an approximation to the physics of transport by drift and orbit scattering effects, which are still under development and called 'advanced neoclassical' theory. See the reference under NCLASS equations for a more complete discussion of the model.
The signs of the magnetic fields and flows are consistent with the right-handed cylindrical coordinate system (R,f,Z) where R is measured from the major axis and Z is up. So the toroidal angle f runs counter-clockwise when viewed from the top of the torus. The same toroidal angle f is used for the right-handed (r,q,f) curvilinear coordinates so q runs downward on the outside midplane.
Most TFTR and DIII-D data have negative toroidal magnetic fields and positive toroidal currents and poloidal magnetic fields. The safety factor 'q' should be negative (as in the test case) if it is to be consistent with this definition of coordinates, but it is usually reported as positive in databases by long-standing convention. JET with its iron core typically operates with toroidal field and current in the same direction. NCLASS gives correct signs for the flows and radial electric field if the input obeys the above sign conventions.
a0
.
The grid units are represented as 'rho' inside NCLASS.
To resolve the axial singularity, choose r proportional to the poloidal flux, Y. See the discussion under NCLASS_PT cases.
The form of the diffusion term in the continuity equation is assumed to not contain <grad(r)2> or similar geometric factors. This is because the geometry is built into the calculation of the fluxes and transport coefficients through the drifts, which are related to the required flux surface averaged input quantities. If the user has <grad(r)2> in the transport equations, then the transport coefficients from NCLASS must be divided by that same factor.
No diffusivities or convective velocities are presently returned
for the external force (component 5). The user can create them
from the returned flux of species s,
gfl_s(5,s)
.
The addition of a radial electric field, -F', does not affect the radial particle or conduction fluxes. The forces from F' are included in NCLASS. Their lack of influence on the poloidal rotation velocities and fluxes can be used as a check on the computations. A radial electric field gradient does, however, modify the flows and fluxes through 'orbit squeezing.'
For the continuity equation,
dni/dt = 1/V'[V'Gi]' + Si, i = species, ' = d/drthe radial fluxes are returned from NCLASS in three different forms:
gfl_s
, broken into five components:
1) banana-plateau flux driven by p' and T' forces,
2) Pfirsch-Schlüter flux,
3) classical flux,
4) banana-plateau flux driven by the parallel electric
field (<E.B>) force, and
5) banana-plateau flux driven by the external force:
Gi = Sk gfl_s(k,i), k = 1,5 components
dn_s
, plus three
contributions to the convective velocity:
1) veb_s
driven by the <E.B> force
(neoclassical pinch),
2) vn_s
driven by all other off-diagonal
p' and T' forces, and
3) the velocity, calculated externally from the flux
gfl_s(5,i)
, driven by the external force,
fex_iz
:
Gi/ni = - dn_s(i) ni'/ni + vn_s(i) + veb_s(i) + gfl_s(5,i)/ni
dt_ss
, and pressure gradients,
dp_ss
, plus the <E.B> and external force components:
Gi/ni = Sj [- dt_ss(j,i) Tj'/Tj - dp_ss(j,i) pj'/pj] + veb_s(i) + gfl_s(5,i)/ni
The same comment about including a geometric factor in the particle fluxes also applies to the conduction fluxes.
If all the ions have the the same temperature and the user wants to combine the conduction heat fluxes into a single effective ion thermal conductivity plus a velocity component (as in the second option above for the particle fluxes), that can be done externally to NCLASS. This reconstruction is done in the stand-alone driver, NCLASS_PT, and presented as the second form below).
The energy balance equation
1.5 d(niTi)/dt = 1/V'[V'(Qi + 2.5 TiGi)]' + ...The three forms of the conduction heat fluxes are:
qfl_s
, for each species broken
into the same components as the particle flux:
Qi = Sk qfl_s(k,i), k=1,5 components
dq_s
, as
well as the corresponding convective component, vq_s
.
In the equations below, subscript a
represents an
'average' ion, and is obtained by summing over all ion species.
Qa
can be used to evaluate
va
, as is done in NCLASS_PT. Below,
electrons are assumed to be the first species, although in both
NCLASS_PT and NCLASS isotopes can be in any order.
Qe = -neceTe' + neTeve ce = chit_ss(1,1) + chip_ss(1,1) ve = Qe/neTe + ceTe'/Te Qa = -nacaTi' + naTava ca = Si ni[chit_ss(i,i) + chip_ss(i,i)]/na va = Qa/naTi + caTi'/Ti na = Si ni Qa = Si Qi
Qi = niTi{Sj [-chit_ss(j,i) Tj'/Tj - chip_ss(j,i) pj'/pj] + qeb_s(i) + qfl_s(5,i)/(niTi)}
p_bsjb
. The bootstrap current is also
returned as a set of coefficients of the pressure and temperature
gradients of each species, bsjbt_s
and
bsjbp_s
:
<Jbs.B> = p_bsjb <Jbs.B> = Si [bsjbt_s(i)Ti'/Ti + bsjbp_s(i)pi'/pi]
p_exjb
. An example of
such a current is the electron response current from neutral
beam injection, which can be entered through
fex_iz
.
<Jex.B> = p_exjb
NCLASS_
extension - NCLASS module and driver
RARRAY_
extension - array manipulation routines
WRITE_
extension - output routines
U_
extension - utility routines
mx_
extension - dimensioning parameter
c_
extension - constant (profile-independent)
quantity
k_
extension - option
p_
extension - single point (profile dependent)
variable
z_
extension - physical or conversion constant
_iz
- charge state array
_s
or
descriptor_ss
- singly or doubly dimensioned
species array, respectively
_i
or
descriptor_ii
- singly or doubly dimensioned
isotope array, respectively
k_order-order of v moments to be solved (-) =2 u and q =3 u, q, and u2 =else error k_potato-option to include potato orbits (-) =0 off =else on
m_i-number of isotopes (1 < m_i < mx_mi+1) m_z-highest charge state of all species (0 < m_z < mx_mz+1)
c_den-density cutoff below which species is ignored (/m3) c_potb-kappa(0)*Bt(0)/[2*q(0)2] (T) c_potl-q(0)*R(0) (m)
p_b2-<B2> (T2) p_bm2-<1/B2> (/T2) p_eb-<E.B> (V*T/m) p_fhat-mu_0*F/(dPsi/dr) (rho/m) p_fm(3)-poloidal moments of geometric factor for PS viscosity (-) p_ft-trapped fraction (-) p_grbm2-<grad(rho)2/B2> (rho2/m2/T2) p_grphi-radial electric field Phi' (V/rho) p_gr2phi-radial electric field gradient Psi'(Phi'/Psi')' (V/rho2) p_ngrth-<n.grad(Theta)> (1/m)
amu_i(i)-atomic mass number of i (-) grt_i(i)-temperature gradient of i (keV/rho) temp_i(i)-temperature of i (keV)
den_iz(i,z)-density of i,z (/m3) fex_iz(3,i,z)-moments of external parallel force on i,z (T*j/m3) grp_iz(i,z)-pressure gradient of i,z (keV/m3/rho)
iflag-warning and error flag =-4 warning: no viscosity =-3 warning: no banana viscosity =-2 warning: no Pfirsch-Schluter viscosity =-1 warning: no potato orbit viscosity =0 no warnings or errors =1 error: order of v moments to be solved must be 2 or 3 =2 error: number of species must be 1< m_i < mx_mi+1 =3 error: number of species must be 0 < m_z < mx_mz+1 =4 error: number of species must be 1 < m_s < mx_ms+1 =5 error: inversion of flow matrix failed =6 error: trapped fraction must be 0.0.le.p_ft.le.1.0
m_s-number of species (1 < m_s < mx_ms+1)
p_bsjb-<J_bs.B> (A*T/m2) p_etap-parallel electrical resistivity (Ohm*m) p_exjb-<J_ex.B> current response to fex_iz (A*T/m2)
calm_i(3,3,i)-tp eff friction matrix for i (kg/m3/s) caln_ii(3,3,i1,i2)-fp eff friction matrix for i1 on i2 (kg/m3/s) capm_ii(3,3,i1,i2)-test part (tp) friction matrix for i1 on i2 (-) capn_ii(3,3,i1,i2)-field part (fp) friction matrix for i1 on i2 (-)
bsjbp_s(s)-<J_bs.B> driven by unit p'/p of s (A*T*rho/m3) bsjbt_s(s)-<J_bs.B> driven by unit T'/T of s (A*T*rho/m3) chip_ss(s1,s2)-heat cond coefficient of s2 on p'/p of s1 (rho2/s) chit_ss(s1,s2)-heat cond coefficient of s2 on T'/T of s1 (rho2/s) dn_s(s)-diffusion coefficient (diag comp) of s (rho2/s) dp_ss(s1,s2)-diffusion coefficient of s2 on p'/p of s1 (rho2/s) dt_ss(s1,s2)-diffusion coefficient of s2 on T'/T of s1 (rho2/s) gfl_s(m,s)-radial particle flux comps of s (rho/m3/s) m=1, banana-plateau, p' and T' m=2, Pfirsch-Schlüter m=3, classical m=4, banana-plateau, <E.B> m=5, banana-plateau, external parallel force fex_iz jm_s(s)-isotope number of s (-) jz_s(s)-charge state of s (-) qeb_s(s)-<E.B> heat convection velocity of s (rho/s) qfl_s(m,s)-radial heat conduction flux comps of s (W*rho/m3) m=1, banana-plateau, p' and T' m=2, Pfirsch-Schlüter m=3, classical m=4, banana-plateau, <E.B> m=5, banana-plateau, external parallel force fex_iz sqz_s(s)-orbit squeezing factor for s (-) upar_s(3,m,s)-parallel flow of s from force m (T*m/s) m=1, p' and T' m=2, <E.B> m=3, fex_iz utheta_s(3,m,s)-poloidal flow of s from force m (m/s/T) m=1, p', T' m=2, <E.B> m=3, fex_iz vn_s(s)-convection velocity (off diag comps-p', T') of s (rho/s) veb_s(s)-<E.B> particle convection velocity of s (rho/s) xi_s(s)-charge weighted density factor of s (-) ymu_s(s)-normalized viscosity for s (kg/m3/s)
amnt_ii(i1,i2)-eff relaxation rate for i1 on i2 (kg/m3/s) vt_i(i)-thermal velocity of i (m/s)
tau_ss(s1,s2)-90 degree scattering time of s1 on s2 (s) ykb_s(s)-banana viscosity for s (kg/m3/s) ykp_s(s)-Pfirsch-Schlüter viscosity for s (kg/m3/s) ykpo_s(s)-potato viscosity for s (kg/m3/s) ykpop_s(s)-potato-plateau viscosity for s (kg/m3/s) ynud_s(s)-pitch angle diffusion rate for s (/s) ynut_s(s)-anisotropy relaxation rate for s (/s) ynuti_s(3,s)-PS anisotropy relaxation rates for s (/s)
k_out-option for output to nout (-) =1 errors only =2 errors and results =else no output
p_eps-inverse aspect ratio (-) p_q-safety factor (-)
r0-major radius (m) a0-minor radius, scale length for rho (m) e0-axial elongation (-) bt0-axial toroidal field (T) q0-axial safety factor (-)
Parameters FORTRAN Source NCLASS_PT_DR mx_mi,mx_ms,mx_mz NCLASS_PT_DR.FOR | |--NCLASS NCLASS_MOD.FOR | | | |--NCLASS_MN mx_mi NCLASS_MOD.FOR | | | |--NCLASS_TAU mx_mi,mx_ms,mx_mz NCLASS_MOD.FOR | | | | | |--RARRAY_ZERO RARRAY_ZERO.FOR | | | |--NCLASS_MU mx_mi,mx_ms,mx_mz NCLASS_MOD.FOR | | | | | |--NCLASS_K mx_mi,mx_ms NCLASS_MOD.FOR | | | | | | | |--NCLASS_NU mx_mi,mx_ms NCLASS_MOD.FOR | | | | | | | | | |--U_ERF U_ERF.FOR | | | | | | | | | |--RARRAY_ZERO RARRAY_ZERO.FOR | | | | | | | |--RARRAY_ZERO RARRAY_ZERO.FOR | | | | | |--RARRAY_ZERO RARRAY_ZERO.FOR | | | |--NCLASS_FLOW mx_mi,mx_ms,mx_mz NCLASS_MOD.FOR | | | |--U_LU_DECOMP U_LU_DECOMP.FOR | | | |--U_LU_BACKSUB U_LU_BACKSUB.FOR | | | |--RARRAY_ZERO RARRAY_ZERO.FOR | |--RARRAY_COPY RARRAY_COPY.FOR | |--RARRAY_SUM RARRAY_SUM.FOR | |--RARRAY_ZERO RARRAY_ZERO.FOR | |--WRITE_IR WRITE_MOD.FOR | |--WRITE_LINE WRITE_MOD.FOR | |--WRITE_LINE_IR WRITE_MOD.FOR