PB3D flowchart

In figure 1, the top-level flowchart of PB3D is shown (adapted from [16]):

Figure 1 : basic PB3D flowchart

Roughly speaking, the main body of the code is taken up by four driver routines:

input driver
equilibrium driver
perturbation driver
solution driver

The input driver sets the other drivers up by reading input options as well as quantities from equilibrium codes, and converting them to the internal PB3D format.

The equilibrium driver uses these input variables to set up equilibrium quantities.

On the one hand these consist of flux equilibrium quantities (stored in a eq_vars.eq_1_type) that only depend on the flux coordinates, such as the pressure \(p\), the safety factor \(q\), etc.
On the other hand, they consist of metric equilibrium quantities (stored in a eq_vars.eq_2_type) that depend on the angular position as well, such as the metric quantities of the grid \(g_{ij} = \vec{e}_i \cdot \vec{e}_j\) and \(h_{ij} = \nabla u^i \cdot \nabla u^j\), as well as derived quantities such as the parallel current \(\sigma = \frac{\vec{B}\cdot\vec{J}}{B^2}\), etc. [15]

The perturbation driver uses both types of equation variables to calculate perturbation quantites. Again, there are two types.

The first type of variables are the vectorial perturbation quantities (stored in a x_vars.x_1_type), which contains the geodesic component of the perturbation \(U = \frac{\nabla \times \vec{B}}{\left|\nabla \psi\right|^2} \cdot \vec{\xi}\) where \(\vec{\xi}\) is the perturbation. Also the parallel derivatives of these quantities are calculated.
The other type re the tensorial perturbation quantities (stored in a x_vars.x_2_type), and it combines the vectorial perturbation quantities with the equilibrium quantities to form the tensorial quantities \(\widetilde{PV}_{k,m}^i\) and \(\widetilde{KV}_{k,m}^i\) that represent the perturbed potential energy as well as the kinetic energy of the perturbation [15].

These are also complemented by \(\delta_{k,m}^\text{vac}\), the vacuum contribution to the perturbed potential energy (1).

After integrating them along the magnetic field lines, these magnetic integrals form the building blocks of the discretized matrices \(\overline{\text{A}}\) and \(\overline{\text{B}}\) in the generalized system of eigenvalue equations \(\overline{\text{A}} \vec{X} = \lambda \overline{\text{B}}\vec{X}\).

These matrices are set up in the solution driver, which currently works with the SLEPc library that is build on the PETSc library for sparse linear algebra. After solution, the eigenvalues and -vectors are stored appropriately.

Finally, it is checked whether another Richardson level should be attempted or not in Richardson converged?, and if so, the number of parallel grid points is doubled.

Equilibrium & Perturbation Jobs

The general code structure explained in PB3D flowchart basically comes down to setting up the matrices to be solved for the eigenvalues and eigenvectors. In general, the dimensions of these variables are too large (see the grids in Detailed PB3D Flowchart) for them to fit in memory. A structure of two nested loops is therefore employed in PB3D, called jobs, to keep the memory requirements under a user-defined maximum (see max_tot_mem in Input variables).

Perturbation jobs

perturbation jobs are used to divided the calculation of the tensorial perturbation quantities, which come in pairs of mode numbers. These combinations are blocked, where each block fits in memory. These blocks are the perturbation jobs.

They form the inner loop.

See also: See divide_x_jobs().

Equilibrium jobs

As a result the memory still needs to hold all the relevant equilibrium and vectorial perturbation quantities necessary for the individual perturbation jobs during the whole time these jobs are being calculated. equilibrium jobs are therefore employed in an outer loop to divide these equilibrium and vectorial perturbation quantities in manageable chunks.

To be more precise, the parallel range over which variables are to be integrated in the magnetic integrals (see Detailed PB3D Flowchart), is cut in pieces. The advantage here is that only the results from the sub-integrals have to be saved and updated, not the function values themselves.

See also: See divide_eq_jobs().

Detailed PB3D Flowchart

The different blocks if figure 1 are expanded in table 1. The cells that are typeset in light blue are only performed for the first equilibrium job of either the restart Richardson level, or the first Richardson level if this is not provided.

Table 1. detailed PB3D flowchart
variable type	VMEC (`eq_style = 1`)	HELENA (`eq_style = 2`)
PB3D initialize	the PB3D structure is set up, such as the MPI communicator, etc.
input driver	The arguments are parsed (see parse_args())
	The input files are opened (see open_input() and Input variables)
	The user inputs are processed (see read_input_opts())
	The equilibrium code output is read and converted to the internal PB3D format (see read_input_eq())
	Normalization constants are calculated (see calc_normalization_const()) and the input is normalized (see normalize_input())
	Output files are opened (see open_output() and Code outputs)
	Input variables are written to HDF5 output (see print_output_in())
	Some other variables are broadcasted directly to all the processes using MPI (see broadcast_input_opts())
equilibrium driver	normal part of equilibrium grid is set up, as for the angular part, the field-lines are necessary which are not yet calculated	the equilibrium grid is set up equal to the HELENA output grid
	flux equilibrium quantities (see eq_vars.eq_1_type) are calculated (see eq_ops.calc_eq())
	flux equilibrium quantities are written to HDF5 output (see eq_ops.print_output_eq())
	equilibrium grid is completed by tracing the field lines (see calc_ang_grid_eq_b())	an additional, field-aligned equilibrium grid is set up by tracing the field lines (see setup_grid_eq_b()), and it is written to HDF5 output (see print_output_grid())
	the equilibrium grid is written to output	the equilibrium grid is written to the output
	the metric equilibrium variables (see eq_vars.eq_2_type) are calculated (see eq_ops.calc_eq())	the metric equilibrium variables (see eq_vars.eq_2_type) are calculated (see eq_ops.calc_eq())
		the metric equilibrium variables are written to the output (see eq_ops.print_output_eq())
	store equilibrium quantities necessary for vacuum calculation (see store_vac())	store equilibrium quantities necessary for vacuum calculation (see store_vac())
	possible plot of magnetic field, current density and/or curvature in (field-aligned) equilibrium grid (see b_plot(), j_plot(), kappa_plot())	possible plot of magnetic field, current density and/or curvature in (HELENA) equilibrium grid (see b_plot(), j_plot(), kappa_plot())
	redistribute the equilibrium grid (see redistribute_output_grid())	redistribute the equilibrium grid (see redistribute_output_grid())
	redistribute the flux equilibrium quantities (see eq_ops.redistribute_output_eq())
	redistribute the metric equilibrium quantities (see eq_ops.redistribute_output_eq())	redistribute the metric equilibrium quantities (see eq_ops.redistribute_output_eq())
		The field-aligned equilibrium grid is redistributed (see redistribute_output_grid())
	possible plot of magnetic field lines in flux surfaces (see magn_grid_plot()) if last equilibrium job
perturbation driver		the metric equilibrium quantities are inerpolated in the field-aligned equilibrium grid (see interp_hel_on_grid())
	determine perturbation grid (see setup_grid_x()) and write it to HDF5 ouput (see print_output_grid())	determine perturbation grid (see setup_grid_x()) and write it to HDF5 ouput (see print_output_grid())
		set up field-aligned perturbation grid and write to output
	set up mode number information (see setup_nm_x()) and check them (see check_x_modes()) if first perturbation job
	possibly plot resonance (see resonance_plot()) if first perturbation job
	the vectorial perturbation variables (see x_vars.x_1_type) are calculated (see x_ops.calc_x())	the vectorial perturbation variables (see x_vars.x_1_type) are calculated (see x_ops.calc_x())
		write vectorial perturbation variables to output (see x_ops.print_output_x())
		interpolate vectorial perturbation variables on field-aligned grid (see interp_hel_on_grid())
	calculate tensorial perturbation variables (see x_ops.calc_x())
	calculate magnetic integrals (see calc_magn_ints())
	write magnetic integrals to output (see x_ops.print_output_x())
solution driver	calculate vacuum (see calc_vac())	calculate vacuum (see calc_vac())
	set up the solution grid (see setup_grid_sol()) and write to HDF5 output (see print_output_grid())
	solve system of equations (see solve_ev_system_slepc())
Richardson converged?	check for convergence of Richardson levels by comparing solution eigenvalues (see check_conv())
PB3D finalize	clean up

POST Overview

The POST program for post-processing of PB3D results consists of just one driver.

In the first part, everything is set up (see init_post()):
- This includes initializing the variables written to HDF5 in PB3D.
- Furthermore, the output grids are created
- Also, 1-D flux plots are produced, which are part of the plots that PB3D is able to produce (see plot_flux_q, plot_resonance and plot_magn_grid in table 3 in Code outputs) can be created here.
After that, new outputs are calculated and plot (see run_driver_post()):
- Quantities are calculated on the output grids.
- Using these, plots are first produced that do not depend on the solution vector from PB3D, and are also part of the plots that PB3D is able to produce (see plot_B, plot_J, plot_kappa in table 3 in Code outputs).
- Finally, other plots are produced that do depend on these (see plot_sol_xi, plot_sol_Q, plot_E_rec in table 3 in Code outputs).

POST also knows the concept of equilibrium jobs, but not that of perturbation jobs (see Equilibrium & Perturbation Jobs). This is necessary for larger outputs that quickly overflow memory. The HDF5 output files are updated at every equilibrium job, but, with reasonable safety, they can already be looked at in a visualizer.

Note

This is not yet implemented.