There is some fault detection built into the routines. These distinguish between minor faults, that trigger a warning but do not stop execution, and major faults, that trigger an error that stop execution.
See also
See writo().
The data type real refers to double precision real, defined by num_vars.dp. For integers, standard precision is used.
There are hard-coded min_tol and max_tol in read_input_opts(), which are used to limit the tolerances tol_rich, tol_zero and tol_SLEPC. This could be changed if needed.
If eq_style = 2 (i.e. HELENA), the total maximum memory available must be enough to perform the first Richardson level in a single equilibrium job (see General code structure).

PB3D Inputs

Input file

input parameter explanation default value data type note

concerning solution

n_r_sol number of points in solution grid $100$ int

min_r_sol minimum normalized flux of computational domain $\left(0\ldots 1\right)$ $0.1$ real (1)

max_r_sol maximum normalized flux of computational domain $\left(0\ldots 1\right)$ $1.0$ real (1)

tol_norm tolerance for normal range $0.05$ real (2)

EV_style

style for EV calculation:

$1$	SLEPC solver

$1$

int

concerning field line

alpha_style

style for alpha:

$1$	one field line, many turns
$2$	many field lines, one turn

$2$ (VMEC),

$1$ (HELENA)

int

n_alpha number of field lines for alpha_style $2$ $10$ (VMEC) int

alpha field line label $\left[\pi\right]$ for alpha_style $1$ $0$ real

min_alpha minimum field line label $\alpha$ $\left[\pi\right]$ for alpha_style $2$ $0$ int

max_alpha

maximum field line label

$\alpha$

$\left[\pi\right]$ for alpha_style

$2$

int

concerning perturbation

min_n_par_X minimum number of parallel points for integration $20$ int (3)

min_par_X minimum parallel angle $\left[\pi\right]$ for integration $-4$ int

max_par_X maximum parallel angle $\left[\pi\right]$ for integration $4$ int

prim_X primary mode number in geodesic direction $20$ int (4)

min_sec_X minimum secondary mode number in parallel direction int (4, 5)

max_sec_X maximum secondary mode number in parallel direction int (4, 5)

n_mod_X number of secondary modes in parallel direction $20$ int (4, 5)

use_pol_flux_F

on which flux to base the normal coordinate.

`.true.`	rescaled poloidal flux $\frac{\Psi_\text{pol}}{2 \pi}$
`.false.`	rescaled toroidal flux $\frac{\Psi_\text{tor}}{2 \pi}$

.true.

log

(35)

concerning normalization

rho_0 Normalization factor for density $\rho$ $\frac{kg}{m^3}$ $10^{-7}$ real

R_0 Normalization factor for length scale $R$ $m$ real (6)

pres_0 Normalization factor for pressure $p$ $Pa$ real (6)

psi_0 Normalization factor for flux $\Psi$ $T m^2$ real (6)

B_0 Normalization factor for magnetic field $\vec{B}$ $T$ real (6)

T_0

Normalization factor for time

$s$

real

(6)

concerning input and output

n_sol_requested number of solutions requested $1$ int

retain_all_sol

retain all solutions

`.true.`	unphysical eigenvalues are also retained
`.false.`	unphysical eigenvalues are left out

$2$

log

(7)

plot_resonance plot the safety factor (poloidal flux) or rotational transform (toroidal flux) and the resonant values .false. log (13)

plot_magn_grid plot magnetic field lines in the flux surfaces in 3-D geometry as an animation .false. log (14)

plot_B plot magnetic field $\vec{B} = \nabla \alpha \times \nabla \psi$ $\left[T\right]$ .false. log (15)

plot_J plot current $\vec{J} = \mu_0^{-1} \nabla \times \left(\nabla \alpha \times \nabla \psi\right)$ $\left[\frac{A}{m^2}\right]$ .false. log (16)

plot_kappa plot curvature $\vec{\kappa} = \frac{\vec{B}}{B} \cdot \nabla \frac{\vec{B}}{B}$ $\left[\frac{1}{m}\right]$ and its components $\kappa_n = \frac{\nabla \psi}{\left|\nabla \psi\right|^2} \cdot \vec{\kappa}$ $\left[\frac{1}{T m^2}\right]$ and $\kappa_g = \frac{\nabla \psi \times \vec{B}}{B^2} \cdot \vec{\kappa}$ $\left[\phantom{\cdot}\right]$ .false. log (17)

plot_flux_q plot flux quantities $q$ , $\iota$ , $p$ $\left[Pa\right]$ , $\Psi_\text{pol}$ $\left[T m^2\right]$ and $\Psi_\text{tor}$ $\left[T m^2\right]$ .false. log (18)

n_theta_plot number of points in plots in $\theta$ direction int (19)

n_zeta_plot number of points in plots in $\zeta$ direction int (19)

min_theta_plot minimum of $\theta$ range for plots $1$ int

max_theta_plot minimum of $\theta$ range for plots $3$ int

min_zeta_plot minimum of $\zeta$ range for plots int (20)

max_zeta_plot minimum of $\zeta$ range for plots int (20)

plot_size size of 2-D plots in inch by inch $\left[10,5\right]$ int(2)

ex_plot_style

style how external plots are made

$1$	GnuPlot
$2$	Bokeh for 2-D and Mayavi for 3-D

$1$

int

(21)

concerning Richardson extrapolation

max_it_rich maximum number of Richardson extrapolations $1$ int

tol_rich tolerance for convergence of Richardson extrapolations $1^{-4}$ real

rich_restart_lvl

Richardson level at which to restart a previous simulation

$1$

int

(8)

concerning finding zeros

max_it_zero maximum number of iterations in num_ops.calc_zero_hh() and calc_zero_zhang() $100$ int (9)

tol_zero tolerance for zero finding $1^{-10}$ real (9)

max_nr_backtracks_HH maximum number of times the correction can be relaxed by half in a backtrack for the Householder methods $20$ int (9)

concerning runtime

use_normalization use normalization for the internal calculations .true. log (6)

max_tot_mem maximum total memory available to the simulation $\left[kB\right]$ 6000 real

sol_n_procs Number of processes available for SLEPC solution $(1\ldots N)$ where $N$ is the number of MPI processes N int

norm_disc_prec_eq normal discretization precision for equilibrium quantities $3$ int

norm_disc_prec_X normal discretization precision for perturbation quantities $3$ int

norm_disc_prec_sol normal discretization precision for solution quantities $3$ int

norm_disc_style_sol

style how solution quantities are discretized

$1$	central finite differences
$2$	left finite differences

$2$

int

magn_int_style

style how magnetic integrals are calculated

$1$	trapezoidal rule
$2$	Simpson 3/8 rule

$1$

int

X_grid_style

style how perturbation grid is set up

$1$	identical to equilibrium grid
$2$	identical to solution grid
$3$	enriched version of equilibrium grid (see `max_njq_change`)

$1$ (X_style 1),

$2$ (X_style 2)

int

(5, 36)

max_njq_change maximum change of resonant mode numbers for X_grid_style $3$ $0.49$ real (36)

max_it_SLEPC maximum number of iterations in SLEPC calculation to find eigenvalue $1000$ int

EV_BC eigenvalue used in boundary condition of style BC_style = 1 $1$ real (10)

EV_guess target guess for eigenvalue $-0.3$ real

tol_SLEPC tolerance used in SLEPC for different Richardson levels $[1^{-5},\ldots]$ real(max_it_rich)

rho_style

style how the density

$\rho$ is calculated

$1$	constant `rho_0`

$1$

int

U_style

style how to calculate geodesical perturbation

$U = \vec{\xi} \cdot \frac{\nabla \psi \times \vec{B}}{\left|\nabla \psi\right|^2}$ is calculated

$1$	up to order $\sim \frac{1}{n}$
$2$	up to order $\sim \left(\frac{1}{n}\right)^2$
$3$	up to order $\sim \left(\frac{1}{n}\right)^3$

$3$

int

K_style

style how to calculate kinetic energy

$K = \int \left|\vec{\xi}\right|^2 d \text{V}$

$1$	normalization of full perpendicular component
$2$	normalization of only normal component

$1$

int

BC_style

style how boundary conditions are applied for left and right boundary

$1$	set to zero
$2$	minimization of surface energy through asymmetric finite differences
$3$	minization through extension of grid
$4$	explicit introduction of the surface energy minimization

$[1,2]$

int(2)

(10)

norm_style

style how to calculate normalization

$1$	MISHKA normalization with magnetic field on axis
$2$	COBRA normalization with pressure on axis

$1$

int

(11)

solver_SLEPC_style

style used in SLEPC to solve the eigenvalue equation

$1$	Krylov-Schur with shift-invert
$2$	generalized Davidson with preconditioner

$1$

int

matrix_SLEPC_style

style used for SLEPC matrices

$1$	sparse matrices
$2$	shell matrices

$1$

int

(12)

Command-line inputs

Table 2. PB3D command-line options
input parameter	explanation	note
`test`	test mode Note Debug version only	(28)
`no_plots`	do not produce any output plots
`no_outputs`	do not produce text output	(29)
`do_execute_command_line`	run the external command line shell commands during execution	(30)
`mem_info`	produce memory profiles for each process Note Debug version only	(31)
`no_guess`	do not use a guess
`jump_to_sol`	jump straight to the solution for the Richardson level at which the simulationes are (re)started, possibly for different solution grid parameters such as `n_r_sol`	(33)
`export_HEL`	create a VMEC input file from HELENA variables, possibly adding a toroidal ripple	(34)
`plot_VMEC_modes`	plot decay of VMEC modes
`invert_top_bottom_H`	invert top and bottom of HELENA equilibrium Note Debug version only

POST Inputs

Input file

input parameter explanation default value data type note

concerning solution

min_r_plot minimum normalized flux of output domain $\left(0\ldots 1\right)$ real (24)

max_r_plot

maximum normalized flux of output domain

$\left(0\ldots 1\right)$

real

(24)

concerning input and output

n_sol_plotted number of solutions plotted n_sol_requested int

pert_mult_factor_POST factor by which the grid can be perturbed, with respect to the maximum perturbation of each solution independently $0$ real (25)

plot_resonance plot the safety factor (poloidal flux) or rotational transform (toroidal flux) and the resonant values .false. log (13)

plot_magn_grid plot magnetic field lines in the flux surfaces in 3-D geometry as an animation .false. log (14)

plot_B plot magnetic field $\vec{B} = \nabla \alpha \times \nabla \psi$ $\left[T\right]$ .false. log (15)

plot_J plot current $\vec{J} = \mu_0^{-1} \nabla \times \left(\nabla \alpha \times \nabla \psi\right)$ $\left[\frac{A}{m^2}\right]$ .false. log (16)

plot_kappa plot curvature $\vec{\kappa} = \frac{\vec{B}}{B} \cdot \nabla \frac{\vec{B}}{B}$ $\left[\frac{1}{m}\right]$ and its components $\kappa_n = \frac{\nabla \psi}{\left|\nabla \psi\right|^2} \cdot \vec{\kappa}$ $\left[\frac{1}{T m^2}\right]$ and $\kappa_g = \frac{\nabla \psi \times \vec{B}}{B^2} \cdot \vec{\kappa}$ $\left[\phantom{\cdot}\right]$ .false. log (17)

plot_flux_q plot flux quantities $q$ , $\iota$ , $p$ $\left[Pa\right]$ , $\Psi_\text{pol}$ $\left[T m^2\right]$ and $\Psi_\text{tor}$ $\left[T m^2\right]$ .false. log (18)

plot_sol_xi plot normal mode $\vec{\xi}$ $\left[m\right]$ and its components $X = \frac{\nabla \psi}{B^2} \cdot \vec{\xi}$ $\left[\frac{m^2}{T}\right]$ and $U = \frac{\nabla \psi \times \vec{B}}{\left|\nabla \psi\right|^2} \cdot \vec{\xi}$ $\left[\phantom{\cdot}\right]$ .false. log (26)

plot_sol_Q plot magnetic field perturbation due to the normal mode $\vec{Q} = \nabla \times \left(\vec{\xi} \times \vec{B}\right)$ $\left[T\right]$ and its components $Q_n = \frac{\nabla \psi}{B^2} \cdot \vec{Q}$ $\left[m\right]$ and $Q_g = \frac{\nabla \psi \times \vec{B}}{\left|\nabla \psi\right|^2} \cdot \vec{Q}$ $\left[\frac{T}{m}\right]$ .false. log (26)

plot_E_rec plot energy reconstruction .false. log (27)

n_theta_plot number of points in plots in $\theta$ direction int (19)

n_zeta_plot number of points in plots in $\zeta$ direction int (19)

min_theta_plot minimum of $\theta$ range for plots $1$ int

max_theta_plot minimum of $\theta$ range for plots $3$ int

min_zeta_plot minimum of $\zeta$ range for plots int (20)

max_zeta_plot minimum of $\zeta$ range for plots int (20)

plot_size size of 2-D plots in inch by inch $\left[10,5\right]$ int(2)

ex_plot_style

style how external plots are made

$1$	GnuPlot
$2$	Bokeh for 2-D and Mayavi for 3-D

$1$

int

(21)

concerning Richardson extrapolation

PB3D_rich_lvl

Richardson level at which to perform post-processing

int

(22)

concerning finding zeros

max_it_zero maximum number of iterations in num_ops.calc_zero_hh() and calc_zero_zhang() $100$ int (9)

tol_zero tolerance for zero finding $1^{-10}$ real (9)

max_nr_backtracks_HH

maximum number of times the correction can be relaxed by half in a backtrack for the Householder methods

$20$

int

(9)

concerning runtime

max_tot_mem maximum total memory available to the simulation $\left[kB\right]$ 6000 real

POST_style

style how post-processing is done

$1$	extended rectangular grid in $\theta$ and $\zeta$
$2$	field-aligned grid used internally by PB3D

$1$

int

(23)

plot_grid_style

style used for output plotting

$0$	3-D plots
$1$	slab plots: $\theta$ and $\zeta$ are sides of a slab
$2$	slab plots with folding of fundamental interval $0\ldots 2 \pi$
$3$	straight cylinder geometry: $\zeta$ is straightened

$0$

int

(23)

Command-line inputs

Table 4. POST command-line options
input parameter	explanation	argument	note
`test`	test mode Note Debug version only		(28)
`no_plots`	do not produce any output plots
`no_outputs`	do not produce text output		(29)
`do_execute_command_line`	run the external command line shell commands during execution		(30)
`mem_info`	produce memory profiles for each process Note Debug version only		(31)
`swap_angles`	swap $\theta$ and $\zeta$ of output grid in `POST_style = 2`
`compare_tor_pos`	compare $\vec{B}$ , $\vec{J}$ and $\vec{\kappa}$ from `plot_B`, `plot_J` and `plot_kappa` at different toroidal positions, as well as the difference in position	`RZ_0(2)`	(32)

Note

The normalized flux is either the poloidal (use_pol_flux_F = .true.) or toroidal (use_pol_flux_F = .false.), divided by its value at the boundary.
Is used in check_x_modes() whether there exists a range in which each of the modes resonates.
This is the number for the first Richardson level. It is doubled for every next level.
Should satisfy the high- $n$ approximation, currently set to > 5 .
If min_sec_X and max_sec_X is set, the user prescribes the secondary mode number range (X_style = 1). If n_mod_X is set, on the other hand, the user prescribes the size of the range, which is automatically chosen (X_style = 2). They cannot be both chosen.
Normalization constants, with exception of rho_0 are set in the routine calc_normalization_const().
In store_results(), it is decided whether the eigenvalue is physical or not, by considering the ratio between real and imaginary part.
This needs to be between 1 and the maximum level that was obtained in a previous solution, plus 1. For technical HDF5 reasons, if Richardson restart is used, the same input fil parameters should be used for max_tot_mem and the same number of processes should be used. Furthermore, it should not exceed max_it_rich.
For the case of num_ops.calc_zero_hh(), at every iteration it is checked whether the correction approximation is better than the original one. If it is not, the correction is relaxed by a factor half. This is done max_nr_backtracks_HH times.
See set_bc().
See calc_normalization_const().
See solve_ev_system_slepc()
See resonance_plot().
See magn_grid_plot().
See b_plot(). The output is done in various coordinate system, as explained in calc_vec_comp().
See j_plot(). The output is done in various coordinate system, as explained in calc_vec_comp().
See kappa_plot(). The output is done in various coordinate system, as explained in calc_vec_comp().
See flux_q_plot().
For eq_style = 1, the equilibrium is axisymmetric, so by default n_theta_plot = 501 and n_zeta_plot = 1; for eq_style = 2, the equilibrium can be 3-D, so n_theta_plot = 201 and n_zeta_plot = 51
For eq_style = 1, the equilibrium is axisymmetric, so by default min_zeta_plot = max_zeta_plot 0; for eq_style = 2, the equilibrium can be 3-D, so min_zeta_plot = 0 and max_zeta_plot = 2.
This refers to the external plots, which are complimentary to the HDF5 plots. See draw_ex().
If output is done that requires a solution, this needs to be between 1 and the maximum level that was obtained in PB3D. See find_max_rich_lvl(). If no solution is found, some of the outputs that do not depend on this still work.Optionally, this can be forced by choosing PB3D_rich_lvl not positive.
The POST_style determines whether the post-processing is done on an extended, new grid, or on the PB3D internal grid. On top of this, through plot_grid_style the user can set the output grid for the post-processing. Note that some combinations of styles do not make sense.
By default, the plot normal range is chosen equal to the solution normal range from min_r_sol and max_r_sol, but it can be restricted. This is useful, for example, to reduce computation time.
The Cartesian components of the result of the division of pert_mult_factor_POST by the maximum value of the solution normal mode amplitude is added to the equilibrium Cartesian grid components to illustrate the perturbation in 3-D. This is done for each solution independently and sequentially. If it is not positive, it is ignored. Note that plot_sol_xi = .true. is necessary for this to be used.
These quantities are calculated in sol_utilities.calc_xuq(). This is done for a time sequene, currently hard-coded in plot_sol_vec(), depending on whether the solution is stable or unstable.
The energy reconstruction is performed in calc_e() if either one of plot_sol_xi, plot_sol_Q or plot_E_rec is true.
Equivalent use is through -t. If this is used, do_execute_command_line should be set to true as well. See generic_tests().
This is used internally as well, to suppress certain output. See run_driver_x() and run_driver_post().
This can fail on local systems, and generally works badly on computational clusters, where the computational nodes often don't have the external tools (see ex_plot_style) installed. By default, an output file is created for the shell commands, which can be executed afterwards. See Shell commands file.
This is done by calling get_mem_usage() and storing the result at every allocation and deallocation of an eq_1_type, eq_2_type, X_1_type, X_2_type, vac_type or sol_type, as well as for input quantities. Furthermore at every user output with writo(), this is also done, and the output is prepended to the output string. Also, the output is written to a file. See Memory usage file.
This is useful to calculate the ripple in these quantities. To use this, it is necessary to have an output grid of only 3 points, where the poloidal projection of the middle point should be the in the middle between the poloidal projections of the two others. Also, you have to provide through RZ_0 the $R$ and $Z$ value of the origin of the geometrical poloidal angle that is used to calculate distances. See grid_utilities.calc_tor_diff().
Jumping straight to the solution is useful if a simulation already calculated all the necessary quantities to set up the solution matrices $\overline{\text{A}}$ and $\overline{\text{B}}$ in $\overline{\text{A}} \vec{X} = \lambda \overline{\text{B}} \vec{X}$ and the solution needs to be redone, for example because a calculation was attempted with a method that did not converge.
This is a rather complicated procedure with many possibilities that is intended for expert use. See create_vmec_input().
Using the toroidal flux is experimental and untested.
Have a look in calc_norm_range(). X_grid_style $1$ is generally too coarse while X_grid_style $2$ is too fine. X_grid_style $3$ was designed to have the both of best world. It only makes sense for X_style 2 (fast), though, as in the enrichment process the fast secondary mode range is used to limit the maximum amount the quantity $n q$ (when using poloidal flux) or $m \iota$ (when using toroidal flux) to max_njq_change between two normal grid points.

Table of Contents

PB3D Inputs

Input file

Command-line inputs

POST Inputs

Input file

Command-line inputs