Code FLOW user manual

Algorithm settings

The solution algorithm uses a red-black, multi-grid approach. The variables to be set are contained in the namelist input_solver. The following variables can be assigned from outside:

n_min is the number of grid points in each direction for the initial grid.
n is the number of grid points in each direction for the final grid.
min_it is the minimum number of iterations for each grid (usually set to 50)
max_it is the maximum number of iterations for each grid (usually set from a few hundreds to ~2000)
accelerate is a logical variable, controlling whether Chebishev acceleration has to be used in the SOR process
fix_orp is used as a constant over relaxation parameter only if the variable accelerate is set to .false.

Obviously, it must be .

A criterion for convergence is specified in the file mgrid.f90 by the parameter eps, and in general should not be changed.

Boundary Conditions

Boundary Conditions are also controlled in the namelist input_solver.

Two different BC’s can be used in the code. The choice is made assigning the value of the variable bc_type:

bc_type = 1 -> ψ is assigned to be 0 in all point outside the domain of integration
bc_type = 2 -> ψ is interpolated in the grid points immediately outside the computational domain
bc_type = 3 -> ψ is interpolated using a linear algorithm involving four points
bc_type = 2 -> bc_type=1 is used for coarse grids, bc_type=3 for fine grids

If bc_type = 2 there is a numerical stability problem near the boundary. This is addressed assigning a limit value for the over relaxation parameter. The value is assigned in the variable max_orp. If bc_type = 1 the value will not be used. Some numerical tests show that 1.5 is a “safe” value to assume for this variable.

Initial constants

Several constants are specified in the namelist input_constants. As for geometrical constants:

Rmajor is the major radius of the torus
x_size and z_size are the dimensions of the computational grid

Other programming constants:

eq_type specifies the closure to be used. 1 is for MHD, 3 for kinetic with toroidal flow only
eq3_opt is an option for the kinetic closure. If it is equal to 1, temperatures are assigned as free functions of psi, if equal to 2, pressures are assigned. The value is irrelevant if eq_type is not 3.
numerical is a logical constant, specifying whether any function of psi is assigned through a numerical input (i.e., tables). If numerical is set to .false., no function is assigned from a numerical input. Otherwise, some or all functions will be assigned from tables.
broot is a constant specifying the type of equilibrium to be sought. For NSTX equilibria with toroidal flow, this will always be equal to 1.
mass and eV are the ion mass and the elementary charge

Magnetic field data

The magnetic field data are contained in the namelist input_magnetic. is assigned, requiring the following input.

B_f in the vacuum is given by the variable b_phi_zero

The variable F_opt controls the type of equation used for F(ψ).

If F_opt is equal to 1, then

(1)

F_plamsa/F_vacuum is assigned in the variable Fc_o_Fv

k is assigned in the namelist as well.

If F_opt is equal to 2, then

(2)

η_P is assigned in the variable eta_P

If F_opt is equal to 3, an analytic function for F will be used. This must be specified by the user in the function bzero.

mu_mag is the permeability of free space.

A different equation for F can be specified in the function bzero. The derivative dbzerodpsi must be modified as well.

Flow data

The flow data are contained in the namelist input_flow. Values for the sonic Mach numbers can be assigned arbitrarily. If eq_type = 3, the poloidal flow is automatically set equal to 0.

The toroidal flow is defined by:

(3)

The constants M_j^MAX, a_mphi and M_j^min are specified in the namelist.

If a different equation for M_j is required, it must be specified in the function mach_phi. The first derivative must be specified consistently in the function dmach_phidpsi.

If poloidal flow is needed, the functions mach_theta and dmach_thetadpsi must be specified in the source code.

Pressure and density profiles

Pressure and density are specified in the namelist input_p_d_profile. If the values in the namelist input_constants are in metric units, the values in this file will be in metric units as well.

gamma is the ratio of specific heats.

Pressure and density are assigned by the equations:

(4)

D_center is assigned in the corresponding variable, while for D_edge it is assigned:

P_center is assigned through the variable β_center. P_edge is derived from:

The exponents in (4) must be assigned as well.

If the system is anisotropic, several options are present.

If the kinetic closure is used, either Quasi-Pressures or Quasi-Temperatures can be assigned. In the first case (eq3_opt = 2), both parallel and perpendicular pressures are assigned as in (4). Parallel and perpendicular input betas are therefore required. The edge values are obtained from

and .

If eq3_opt is equal to 1, parallel temperature is assigned as in (4), while perpendicular temperature is obtained from:

(5)

where

(6)

and is a measure of the anisotropy of the system.

The relevant variables to be assigned are tparcenter and tpae_o_tpac. theteps gives the function:

The variables tperpcenter and tpee_o_tpec are not used in this version of the code.

Plasma shape and triangularity

Triangularity is controlled in the namelist input_triangularity. The equation for the minor radius r is assigned by the parameter tri_type, as detailed in the following table:

If tri_type=0, a_elps and b_elps need to be specified, and the boundary is given by:

The coefficients k, d and r_m for the options 1 to 3 must be specified in the module triangularity, and can in general be different above and below the mid-plane.

If tri_type is 4, the radius is given in the file r0_spline2.dat as a list of q,r points. The first and last point must be coincident (with angles differing by 2p).

If tri_type is 6, the program will look for a set of Fourier coefficients for R and Z as functions of the angle. These are the same coefficients given in the shot file #.flo. The coefficients (numerical values only) must be stored in the file R_Z.dat. The file must also contain the geometrical major radius. The boundary will be assigned using the Ezspline library.

If tri_type is 7, the same boundary as in tri_type=6 will be assigned, but in this case, R0 must be defined as shown in the attached picture, basically as the R where . This option does not use any external library.

Numerical input

The namelist input_numerical specifies which free functions of ψ are assigned from numerical data, as summarized in the following table:

flag	function	data file
Numerical_n	D(ψ)	n.dat
Numerical_p_iso	P(ψ)	p_iso.dat
Numerical_p_par	P(ψ)	p_par.dat
Numerical_p_perp	P(ψ)	p_perp.dat
Numerical_F	B₀(ψ)	F.dat
numerical_omega	M_φ(ψ)	omega.dat

The files must contain a list of (ψ,function) values, with ψ ranging from 0 to 1, where 0 corresponds to the boundary and 1 to the magnetic axis. The points do not need to be equally spaced.

Again, if the variable numerical in the namelist input_constants is set to .false., no numerical input will be used, regardless of the content of the namelist input_numerical.

Output

While the code is running, the solution process is shown on the screen in the following way:

Figure 7: running output

Once the run is completed, the final value of ψ_center is printed on the screen, together with the position of the magnetic axis. The output of the code consist in files containing the calculated values for eache of the following quantities:

Psi

Rho

P_par

P_perp

Beta

Beta_par

Beta_perp

Mach_Alfven poloidal

Mach_phi (sonic Mach number)

V_phi

V_poloidal

V_r

V_z

B_phi

B_r

B_z

B_pol (poloidal magnetic field)

B_field (total magnetic field)

j_phi (toroidal current)

Residual (solution error)

Temperature

T_par

T_perp

For each quantity, two different files are saved:

1. a tecplot (.plt) file, containing the 2D output of the variable. The data are arranged as (R,Z,variable), with Z varying faster.

2. a .txt file, containing a line cut of the variable along the midplane