|
PB3D
[2.45]
Ideal linear high-n MHD stability in 3-D
|
This page contains a hands-on tutorial to get you used to running PB3D to calculate 3-D ideal linear high-n MHD stability, and POST to post-process the results.
Be sure to check out the Installation instruction first.
If at any point a routine is not clear, don't hesitate to consult the Procedures and Variables page.
PB3D calculates stability of an equilibrium magnetic toroidal configuration. It is written in a modular way to isolate this main task from the actual code that was used to obtain the equilibrium configuration, and the format that this code may use. Currently, two formats are possible:
VMEC is a general 3-D code built on an interesting idea: The poloidal coordinate \(\theta\) is deformed in such as a way to make the Fourier basis that is used to describe the cylindrical coordinates \(R(s,\theta,\varphi)\) and \(Z(s,\theta,\varphi)\) as narrow as possible, where \(s\) is the flux label ( \(s=\sqrt{\Psi_\text{tor}/\Psi_{\text{tor,edge}}})\) and \(\varphi\) is the regular cylindrical angle (1). This is done by introducing a deformation factor \(\lambda(s,\theta,\varphi)\) defined so that the straight-field line coordinate \(\theta_\text{F}\) is given by \(\theta_\text{F} = \theta + \lambda \).
The fact that \(R\), \(Z\) and \(\lambda\) are given as a double Fourier series in \(\theta\) and \(\zeta\) ensures that the angular derivatives can be done analytically. This property is used in PB3D. Apart from that, the flux functions that describe the profiles on the flux surfaces, the rotational transform \(\iota(s)\) and pressure \(p(s)\), are also returned.
For HELENA, the situation is different, as this code considers axisymmetric equilibria only. Its output consists, aside from flux functions safety factor \(q(\psi)\), presure \(p(\psi)\) and toroidal covariant magnetic field \(F(\psi)\), where \(\psi=\Psi_\text{pol}/\Psi_{\text{pol,edge}}\) is the normalized poloidal angle, of a 2-D table for the upper metric factors \(h^{\psi\psi}\) and \(h^{\psi\theta}\) and the lower metric factor \(g_{\zeta\zeta}\). HELENA uses the same angular coordinates \(\theta\) and \(\zeta = -\varphi\) as PB3D (see 1).
To get the appropriate equilibrium code results into PB3D, it suffices to place the output files in the same directory in which PB3D is run:
fort.12),wout file,and to specify their names as the second input argument, the first one being the file with user inputs (see Setting Up the Input File).
wout file has to be in Netcdf's .nc format. The old .txt format is no longer supported.In any case, if you forget how to run the PB3D (or POST) programs, you can always just type
or even just
IAS and B0 to the mapping file, after the variable JS0 and after the variable RAXIS. Some versions do not write these by default. If yours is one of the above, the solution is to change mapping and recompile and re-run HELENA.It is always important that the equilibrium has been properly calculated, well-converged, etc.
It is possible to use PB3D to check this. This will be quickly described in this optional section.
The checks described below are only possible with the code compiled with the debug flag on (see Installation). Using the debug version, the tests can be triggered using the –test option when running PB3D (or POST), and using the interactive user interface.
They generally output HDF5 files that can be read with ParaView or VisIt to check the difference between two plots. After having a look in the code (search for ltest), this should be self-explanatory.
B_F and on the pressure balance consistency are vital.| tests | purpose |
|---|---|
calculation of g_V | Check whether equilibrium metric coefficients calculated for VMEC equilibria through application of transformation matrices match with a direct implementation of the VMEC metric coefficients. Note that the reason why the transformation matrices are used, is that derivatives can be easily calculated as well, with high precision. |
calculation of jac_V | Similar to the test on the calculation of g_V, but this time for the Jacobian of VMEC coordinates. |
| consistency of HELENA metric factors | Check whether the metric factors returned by HELENA are consistent with the ones calculated from the flux quantities, \(R\) and \(Z\). Basically, this is a measure for the accuracy with which the pressure balance is satisfied. |
| harmonic content of HELENA output | Analyze harmonic content of HELENA variables. |
calculation of D1 D2 h_H | Check whether the cross-derivatives are correctly and accurately calculated. |
calculation of T_EF | Check is \(\mathcal{T}_\text{E}^\text{F}\) complies with the (unpublished) theory. |
calculation of jac_F | Similar to the test on the calculation of jac_V, but this time for the Jacobian of Flux coordinates. |
calculation of B_F | Similar to the test on the calculation of jac_F, but this time for the covariant components and magnitude of the magnetic field in Flux coordinates. |
| consistency of pressure balance | Check how well the pressure balance is satisfied. This is a very crucial requirement for the correct calculation of stability: If the pressure balance is not satisfied up to high precision, there is not much for PB3D to do. |
calculation of DU | Check how well the poloidal derivative of the geodesical component of the minimizing perturbation, \(U\), is calculated by comparing it to the numerical derivative. |
For the specific case of VMEC equilibria, finally, there is the additional –plot_VMEC_modes option that produces outputs of the harmonics of the different quantities: \(R\), \(Z\), \(\lambda\), \(\mathcal{J}\) and \(\vec{B}\).
It is very important that these are well-converged, as the internal consistency depends strongly on this. Therefore, if this is not the case, be prepared to receive unaccurate stability results.
Now that we have the equilibrium files, it is time to have a look at the input file.
Have a look:
min_n_par_X is the number of points in the parallel grid.min_par_X to max_par_X. As this is an axisymmetric example, they can simply be chosen to be a fundamental period of \( 2\pi \). For 3-D equilibria, this number has to be increased until convergence is reached [16].alpha indicates the field line at which the perturbations are situated.n_r_sol is the number of points in the solution grid, that ranges from min_r_sol to max_r_sol in the PB3D normal coordinate (see num_vars.use_pol_flux_f).n_mod_X is the number of resonating Harmonics with prim_X the primary mode number.BC_style(1:2) indicates the style with which the boundary conditions are implemented deep in the plasma and on the plasma edge.max_tot_mem is the total virtual memory that PB3D may consume.use_normalization indicates whether normalization is used (recommended), but it is useful to turn this off for debugging.BC_style(2) is currently set to 1, indicating that the mode amplitude is zero at the edge. In the upcoming version of PB3D, the vacuum perturbation on the potential energy will be taken into account, which enables other boundary condition styles at the plasma edge.There are also some extra options listed (but not loaded in PB3D as they com after the end of the inputdata_PB3D block of the Fortran namelist) below:
They include:
rich_restart_lvl is used to restart a simulation at a certain Richardson level. The Richardson level 1 refers to the first simulation, with the number of parallel points used in the calculation of line integrals along the magnetic field set by min_n_par_X (see inputs_PB3D_file_tab). For every subsequent Richardson level, this number is approximately doubled in order to gain a grid size that is twice as fine. The results from these different levels are extrapolated to a theoretically infinitely fine grid. These operations are done in module rich_ops.tol_rich is the tolerance at which Richardson extrapolation is considered to have converged.tol_SLEPC is an array of tolerances used for the different SLEPc eigenvalue solutions.There are many, many more input parameters. A short description of these can be found in Input variables.
This tutorial provides a basic VMEC equilibrium file to experiment with, called cbm18, which represents a simple circular test case designed to be ballooning unstable. [5] Some properties are given in table 1.
| toroidal current \(I_\text{tot}\) | \(1.614MA\) |
| \(\beta_\text{pol}\) | \(0.443 \) |
| aspect ratio \(\epsilon\) | \(1.5 \) |
Put this equilibrium file, as well as the input file in a separate run directory; let's call it cbm18.
In order to produce outputs, PB3D needs three sub-directories in the run directory, called Data/, Scripts/ and Plots/.
Data/, datafiles are stored that are used for plots.Scripts/, scripts are generated that use the datafiles from Data.Plots/, the resulting plots are located, both for external plots, and for HDF5 plots that can be read later in visualization software such as VisIt or ParaView. In Output Plots, this is explained in more detail.If these directories are not present, you will receive error message, so go ahead and create them.
We can now run the code using mpirun. For example, with
4 processes are used. You can find the VMEC file here.
There are optional run-time options that can be triggered, such as –jump_to_sol, which can be used to easily change solution preferences, such as discretization order, etc. For an overview, see Command-line inputs; they are not treated here.
You will see formatted output on the screen, with indentation that consistent indentation according to the depth of a routine in the program. The same output is also written to an output file, unformatted but still indented (see Log file).
The resulting eigenvalues can also be read in a separate output file (see Eigenvalue summary file).
Take a look at these output files and shuffle a bit with the input parameters if you want. For example, change n_mod_X to another value in order to get a different number Fourier harmonics, or n_r_sol to change the number of normal grid points in the solution. The solution should not change a lot if you have enough harmonics and normal grid poins.
min_r_sol to be slightly larger than 0; 0.01 for example.PB3D comes with a bunch of extra run tools, which can greatly improve user experience.
In table 2, an overview is given:
run | Run script that can be used for PB3D and POST with an extensive range of features including
|
extract_jobs_info | Extracts the eigenvalues from the results of a parameter scan produced by run with input variable modifications. |
inspect_SLURM_jobs | Finds the directories of all the jobs currently running in SLURM and in every directory, it gives the tail of PB3D.out. The user can decide to cancel the job. |
POST provides a way to visualize the results from PB3D. Many options are possible, and they can be activated using plot_ flags (see Input file).
An example is given in POST.input
It can be seen that apart from the plot_ flags, there are also
min_r_plot and max_r_plot, which are the direct equivalent of min_r_sol and max_r_sol in PB3D.min_theta_plot, max_theta_plot and n_theta_plot which indicate the range of the poloidal variable that is plot, as well as the the number of point. There is an equivalent for the toroidal variable.POST_style that allows the user to change between different styles, such as whether the output grid is chosen to be fill the torus in a regular way, for example, or whether it follows the PB3D equilibrium grid.plot_grid_style, subsequently, can be used to change the way in which these grids are produced. For example, they could be changed to slab plots or plots with a straightened toroidal angle.pert_mult_factor_POST can be used to deform the output grids by the perturbation itself.Also here, there are also some extra options listed below:
They include:
ex_plot_style to indicate which program is to be used for external plotting. Changing it to 2 would engage Bokeh and Mayavi.PB3D_rich_lvl can be optionally employed to specify explicitely which Richardson level to use for the post-processing.All of these options are explained in Input file. Again, there are also run-time options that can be found in Command-line inputs.
Now, run the POST program, for example using 4 MPI processes:
An interesting output is the energy reconstruction, which is activated with plot_E_rec. This calculation reinserts the solution normal mode into the perturbed plasma potential terms, in order to check whether the Rayleigh quotient still holds, and to see which terms are dominant and which ones are less important. The results are stored in a separate files as well (see Energy reconstruction file).
If –do_execute_command_line was not explicitely used, the actual external plot programs are not used yet, but their usage is given in a shell commands script file (see Shell commands file). Run this shell with
If everything is okay, you should now see output in the Plots/ folder. If an error occurs, please verfiy whether you have appropriate versions of the external plot programs (i.e. GnuPlot, Bokeh, ...)
Also, open some of the .xmf files in VisIt or ParaView to explore these plots.
To finish this tutorial, feel free to adjust some input parameters. For example, you can change the output grid to become 3-D by setting n_zeta_plot and max_zeta_plot, or you could deform it using pert_mult_factor_POST.