|
PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
|
Operations on the equilibrium variables. More...
Interfaces and Types | |
| interface | calc_eq |
| Calculate the equilibrium quantities on a grid determined by straight field lines. More... | |
| interface | print_output_eq |
| Print equilibrium quantities to an output file: More... | |
| interface | redistribute_output_eq |
| Redistribute the equilibrium variables, but only the Flux variables are saved. More... | |
| interface | calc_rzl |
| Calculate \(R\), \(Z\) & \(\lambda\) and derivatives in VMEC coordinates. More... | |
| interface | calc_g_c |
| Calculate the lower metric elements in the C(ylindrical) coordinate system. More... | |
| interface | calc_g_v |
| Calculate the lower metric coefficients in the equilibrium V(MEC) coordinate system. More... | |
| interface | calc_g_h |
| Calculate the lower metric coefficients in the equilibrium H(ELENA) coordinate system. More... | |
| interface | calc_g_f |
| Calculate the metric coefficients in the F(lux) coordinate system. More... | |
| interface | calc_jac_v |
| Calculate \(\mathcal{J}_\text{V}\), the jacobian in V(MEC) coordinates. More... | |
| interface | calc_jac_h |
| Calculate \(\mathcal{J}_\text{H}\), the jacobian in HELENA coordinates. More... | |
| interface | calc_jac_f |
| Calculate \(\mathcal{J}_\text{F}\), the jacobian in Flux coordinates. More... | |
| interface | calc_t_vc |
| Calculate \(\overline{\text{T}}_\text{C}^\text{V}\), the transformation matrix between C(ylindrical) and V(mec) coordinate system. More... | |
| interface | calc_t_vf |
| Calculate \(\overline{\text{T}}_\text{V}^\text{F}\), the transformation matrix between V(MEC) and F(lux) coordinate systems. More... | |
| interface | calc_t_hf |
| Calculate \(\overline{\text{T}}_\text{H}^\text{F}\), the transformation matrix between H(ELENA) and F(lux) coordinate systems. More... | |
Functions/Subroutines | |
| integer function, public | flux_q_plot (grid_eq, eq) |
| Plots the flux quantities in the solution grid. | |
| integer function, public | calc_derived_q (grid_eq, eq_1, eq_2) |
| Calculates derived equilibrium quantities system. | |
| integer function, public | calc_normalization_const () |
| Sets up normalization constants. | |
| subroutine, public | normalize_input () |
| Normalize input quantities. | |
| integer function, public | b_plot (grid_eq, eq_1, eq_2, rich_lvl, plot_fluxes, xyz) |
| Plots the magnetic fields. | |
| integer function, public | j_plot (grid_eq, eq_1, eq_2, rich_lvl, plot_fluxes, xyz) |
| Plots the current. | |
| integer function, public | kappa_plot (grid_eq, eq_1, eq_2, rich_lvl, xyz) |
| Plots the curvature. | |
| integer function, public | delta_r_plot (grid_eq, eq_1, eq_2, xyz, rich_lvl) |
| Plots HALF of the change in the position vectors for 2 different toroidal positions, which can correspond to a ripple. | |
| integer function, public | divide_eq_jobs (n_par_x, arr_size, n_div, n_div_max, n_par_x_base, range_name) |
| Divides the equilibrium jobs. | |
| integer function, public | calc_eq_jobs_lims (n_par_x, n_div) |
Calculate eq_jobs_lims. | |
Variables | |
| logical, public | debug_calc_derived_q = .false. |
| plot debug information for calc_derived_q() | |
| logical, public | debug_j_plot = .false. |
| plot debug information for j_plot() | |
| logical, public | debug_create_vmec_input = .false. |
| plot debug information for create_vmec_input() | |
Operations on the equilibrium variables.
| integer function, public eq_ops::b_plot | ( | type(grid_type), intent(inout) | grid_eq, |
| type(eq_1_type), intent(in) | eq_1, | ||
| type(eq_2_type), intent(in) | eq_2, | ||
| integer, intent(in), optional | rich_lvl, | ||
| logical, intent(in), optional | plot_fluxes, | ||
| real(dp), dimension(:,:,:,:), intent(in), optional | xyz ) |
Plots the magnetic fields.
If multiple equilibrium parallel jobs, every job does its piece, and the results are joined automatically by plot_HDF5.
The outputs are given in contra- and covariant components and magnitude in multiple coordinate systems, as indicated in calc_vec_comp().
The starting point is the fact that the magnetic field is given by
\[\vec{B} = \frac{\vec{e}_{\theta}}{\mathcal{J}}, \]
in F coordinates. The F covariant components are therefore given by
\[B_i = \frac{g_{i3}}{\mathcal{J}} = \frac{\vec{e}_i \cdot \vec{e}_3}{\mathcal{J}}, \]
and the only non-vanishing contravariant component is
\[B^3 = \frac{1}{\mathcal{J}}. \]
These are then all be transformed to the other coordinate systems.
| [in,out] | grid_eq | equilibrium grid |
| [in] | eq_1 | flux equilibrium variables |
| [in] | eq_2 | metric equilibrium variables |
| [in] | rich_lvl | Richardson level |
| [in] | plot_fluxes | plot the fluxes |
| [in] | xyz | X, Y and Z of grid |
Definition at line 5698 of file eq_ops.f90.
| integer function, public eq_ops::calc_derived_q | ( | type(grid_type), intent(in) | grid_eq, |
| type(eq_1_type), intent(in) | eq_1, | ||
| type(eq_2_type), intent(inout), target | eq_2 ) |
Calculates derived equilibrium quantities system.
S kappa_n kappa_g sigma Naive implementations of these quantities, with the exception of S, give numerically very unfavorable results, so special care must be taken in this procedure.
This is an important issue as these derived equilbrium quantities are important building blocks of the perturbed potential energy, including the drivers of instabilities.
The general formulas are derived under the first section, VMEC. For axisymmetric HELENA, however, simplified equations are possible and these are presented afterwards.
The local shear \(S\) is calculated using equation 3.22 from [17] :
\(S = - \frac{1}{\mathcal{J}} \frac{\partial}{\partial \theta} \left(\frac{ h^{\psi\alpha}}{ h^{\psi \psi}}\right)\)
which doesn't pose any particular problems as there are only angular derivatives.
This is calculated using the parallel derivative of the parallel unit vector (i.e. theta if using poloidal flux and zeta if using toroidal flux).
By taking instead the derivative of the covariant basis vector and realizing that the difference with the real curvature lies solely in a component along the parallel direction, which cancels out when taking the normal or geodesic components, the situation becomes easier.
Asuming the poloidal flux is used as normal coordinate, the taks is therefore how to calculate \(\frac{1}{\left|\vec{e}_\theta\right|^2} \frac{\partial}{\partial \theta} \vec{e}_\theta\).
By transforming both the derivative as well as the basis vector to Equilibrium coordinates, this can be written as
\(\sum_{i=2,3} \sum_{j=2,3} \mathcal{T}_\text{F}^\text{E} \left(3,i\right)) \frac{\partial}{\partial u^i_\text{E}} \left( \mathcal{T}_\text{F}^\text{E} \left(3,j\right) \vec{e}_{j,\text{E}}\right)\)
where the summations only run over the angular coordinates because \(\vec{e}_\theta\) never has any component along \(\vec{e}_{\psi,\text{E}}\).
This double summation formula can then be written as a matrix equation
\(\begin{pmatrix}\mathcal{T}_\text{F}^\text{E}\left(3,2\right) & \mathcal{T}_\text{F}^\text{E}\left(3,3\right)\end{pmatrix} \left[ \begin{pmatrix} \frac{\partial}{\partial u^2} \vec{e}_{2} & \frac{\partial}{\partial u^2} \vec{e}_{3} \\ \frac{\partial}{\partial u^3} \vec{e}_{2} & \frac{\partial}{\partial u^3} \vec{e}_{3}\end{pmatrix}_\text{E} \begin{pmatrix}\mathcal{T}_\text{F}^\text{E}\left(3,2\right) \\ \mathcal{T}_\text{F}^\text{E}\left(3,3\right)\end{pmatrix} + \begin{pmatrix} \frac{\partial}{\partial u^2} \mathcal{T}_\text{F}^\text{E}\left(3,2\right) & \frac{\partial}{\partial u^2} \mathcal{T}_\text{F}^\text{E}\left(3,3\right) \\ \frac{\partial}{\partial u^3} \mathcal{T}_\text{F}^\text{E}\left(3,2\right) & \frac{\partial}{\partial u^3} \mathcal{T}_\text{F}^\text{E}\left(3,3\right) \end{pmatrix} \begin{pmatrix}\vec{e}_2 \\ \vec{e}_3\end{pmatrix} \right]\)
which can be represented shorthand as \(\mathcal{T}_\text{F}^\text{E}\left(3,2:3\right) \left[ \mathcal{D}\vec{e}_\text{E} \mathcal{T}_\text{F}^\text{E}\left(3,2:3\right)^T +\mathcal{D} \mathcal{T} \vec{e}_\text{E}^T \right]\)
It is relatively easy to set up the matrix \(\mathcal{D}\vec{e}_\text{E}\) as a function of covariant basis vectors in the Cylindrical coordinate system. The derivatives of the transformation matrix itself can likewise be found.
The steps used in this routine are therefore
An advantage of using this formulation is that no normal derivatives are needed, so that nothing has to implicitely cancel out.
The parallel current is calculated from the shear with the help of equation 3.29 of [17] :
\(\mu_0 \sigma = -\frac{1}{B_\theta} \left[ 2 \frac{\nabla \psi \times \vec{B}}{\left|\nabla \psi\right|^2} \cdot \frac{\partial \left(\nabla \psi\right)}{\partial \theta} + \mathcal{J} \left|\nabla \psi\right|^2 S \right]\) .
where a similar technique can be used as above, for the calculation of the curvature: As
\(\frac{\nabla \psi \times \vec{B}}{\left|\nabla \psi\right|^2} = \frac{B_\theta \vec{e}_\alpha - B_\alpha \vec{e}_\theta} {\mathcal{J} \left|\nabla\psi\right|^2}\) ,
The parallel derivative of \(\nabla \psi\) can be rewritten in terms of contravariant components of derivatives of covariant basis vectors:
\(\left\{ \begin{aligned} \vec{e}_\alpha \cdot \frac{\partial \left(\nabla \psi\right)}{\partial\theta} = - \nabla \psi \cdot \frac{\partial \vec{e}_\theta}{\partial \alpha} \\ \vec{e}_\theta \cdot \frac{\partial \left(\nabla \psi\right)}{\partial\theta} = - \nabla \psi \cdot \frac{\partial \vec{e}_\theta}{\partial \theta} \\ \end{aligned}\right. \) ,
so that the result is:
\(\frac{\nabla \psi \times \vec{B}}{\left|\nabla \psi\right|^2} \cdot \frac{\partial \left(\nabla \psi\right)}{\partial \theta} = \frac{1}{\mathcal{J}^2} \frac{\nabla \psi}{\left|\nabla\psi\right|^2} \cdot \left[ g_{\alpha\theta} \frac{\partial \vec{e}_\theta}{\partial \theta} - g_{\theta\theta} \frac{\partial \vec{e}_\theta}{\partial \alpha} \right] \) ,
which is given by a formula similar to the one used above for the geodesical curvature.
In debug mode, it can be checked whether the current is indeed divergence-free, with the help of equation 3.33 of [17].
\( -2 p' \int_{\theta_0}^\theta \mathcal{J} \kappa_g \text{d}{\theta} = \sigma\left(\theta\right) - \sigma_0\)
and whether a direct, naive implementation of the parallel current from equation 3.26 of [17] agrees with the more accurate results used here:
\(\mu_0 \sigma = \frac{\partial B_\psi}{\partial \alpha} - \frac{\partial B_\alpha}{\partial \psi} - \mu_0 p' \mathcal{J} \frac{B_\alpha}{B_\theta}\) .
The reason why they are generally different is that this implementation relies on the cancellation of large terms.
As \(B_\alpha = F\left(\psi\right)\) and \(\vec{e}_{\alpha,\text{F}} = \vec{e}_{\phi,\text{H}}\), the naive expression for the shear becomes simply
\(\sigma = -F' - \mu_0 p' \frac{F}{B^2}\) ,
which can be used like that.
The calculate the local shear \(S\) is looks like it is best to use equation 3.22 from [17] , just as in the VMEC case:
As a test, however, equation 3.29 of [17] can be used in stead:
\(\mathcal{J}\left|\nabla \psi\right|^2 S + \mu_0 \mathcal{J}B^2 \sigma = - 2 \frac{F}{R} \left( \frac{Z_\theta}{R} + \frac{Z_\theta R_{\theta\theta} - R_\theta Z_{\theta\theta}}{R_\theta^2 + Z_\theta^2} \right)\).
from which \(\sigma\) can be obtained.
Also the curvature expressions can be simplified for axisymmetric situations. The result is given by
\(\kappa_n = \frac{q R}{F} \frac{Z_\theta \left( R_{\theta\theta} - q^2 R\right) - R_\theta Z_{\theta\theta}} {\left(R_\theta^2 + Z_\theta^2 + \left(q R\right)^2\right) \left( R_\theta^2 + Z_\theta^2 \right)} \)
\(\kappa_g = q R \frac{R_\theta \left( 2 \left(R_\theta^2 + Z_\theta^2\right) + \left(qR\right)^2 \right) - R \left(R_\theta R_{\theta\theta} + Z_\theta Z_{\theta\theta}\right)} {\left(R_\theta^2 + Z_\theta^2 + \left(q R\right)^2\right)^2} \)
| [in] | grid_eq | equilibrium grid |
| [in] | eq_1 | flux equilibrium variables |
| [in,out] | eq_2 | metric equilibrium variables |
Definition at line 4286 of file eq_ops.f90.
| integer function, public eq_ops::calc_eq_jobs_lims | ( | integer, intent(in) | n_par_x, |
| integer, intent(in) | n_div ) |
Calculate eq_jobs_lims.
Take into account that every job has to start and end at the start and end of a fundamental integration interval, as discussed in divide_eq_jobs():
magn_int_style=1 (trapezoidal), this is 1 point,magn_int_style=2 (Simpson 3/8), this is 3 points.for POST, there are no Richardson levels, and there has to be overlap of one always, in order to have correct composite integrals of the different regions.
n_par_X passed into this procedure refers to the quantity that is already possibly halved if the Richardson level is higher than 1. This information is then reflected in the eq_jobs_lims, which refer to the local limits, i.e. only the parallel points currently under consideration.| [in] | n_par_x | number of parallel points in this Richardson level |
| [in] | n_div | nr. of divisions of parallel ranges |
Definition at line 6700 of file eq_ops.f90.
| integer function, public eq_ops::calc_normalization_const |
Sets up normalization constants.
eq_style=1):norm_style: eq_style=2):norm_style rho_0 is not given through by the equilibrium codes and should be user-suppliedDefinition at line 5404 of file eq_ops.f90.
| integer function, public eq_ops::delta_r_plot | ( | type(grid_type), intent(inout) | grid_eq, |
| type(eq_1_type), intent(in) | eq_1, | ||
| type(eq_2_type), intent(in) | eq_2, | ||
| real(dp), dimension(:,:,:,:), intent(in) | xyz, | ||
| integer, intent(in), optional | rich_lvl ) |
Plots HALF of the change in the position vectors for 2 different toroidal positions, which can correspond to a ripple.
Also calculates HALF of the relative magnetic perturbation, which also corresponds to a ripple.
Finally, if the output grid contains a fundamental interval \(2\pi\), the proportionality between both is written to a file.
| [in,out] | grid_eq | equilibrium grid |
| [in] | eq_1 | flux equilibrium variables |
| [in] | eq_2 | metric equilibrium variables |
| [in] | xyz | X, Y and Z of grid |
| [in] | rich_lvl | Richardson level |
Definition at line 6171 of file eq_ops.f90.
| integer function, public eq_ops::divide_eq_jobs | ( | integer, intent(in) | n_par_x, |
| integer, dimension(2), intent(in) | arr_size, | ||
| integer, intent(inout) | n_div, | ||
| integer, intent(in), optional | n_div_max, | ||
| integer, intent(in), optional | n_par_x_base, | ||
| character(len=*), intent(in), optional | range_name ) |
Divides the equilibrium jobs.
For PB3D, the entire parallel range has to be calculated, but due to memory limits this has to be split up in pieces. Every piece has to be able to contain the equilibrium variables (see note below), as well as the vectorial perturbation variables. These are later combined into tensorial variables and integrated.
The equilibrium variables have to be operated on to calculate them, which translates to a scale factor mem_scale_fac. However, in the perturbation phase, when they are just used, this scale factor is not needed.
In its most extreme form, the division in equilibrium jobs would be the individual calculation on a fundamental integration integral of the parallel points:
magn_int_style=1 (trapezoidal), this is 1 point,magn_int_style=2 (Simpson 3/8), this is 3 points.For HELENA, the parallel derivatives are calculated discretely, the equilibrium and vectorial perturbation variables are tabulated first in this HELENA grid. This happens in the first Richardson level. In all Richardson levels, afterwards, these variables are interpolated in the angular directions. In this case, therefore, there can be no division of this HELENA output interval for the first Richardson level.
This procedure does the job of dividing the grids setting the global variables eq_jobs_lims.
The integration of the tensorial perturbation variables is adjusted:
In fact, the equilibrium jobs have much in common with the Richardson levels, as is attested by the existence of the routines do_eq() and eq_info(), which are equivalent to do_rich() and rich_info().
In POST, finally, the situation is slightly different for HELENA, as all the requested variables have to fit, including the interpolated variables, as they are stored whereas in PB3D they are not. The parallel range to be taken is then the one of the output grid, including a base range for the variables tabulated on the HELENA grid. Also, for extended output grids, the size of the grid in the secondary angle has to be included in n_par_X (i.e. toroidal when poloidal flux is used and vice versa). Furthermore, multiple equilibrium jobs are allowed.
To this end, optionally, a base number can be provided for n_par_X, that is always added to the number of points in the divided n_par_X.
g_FD, h_FD and jac_FD are counted, as the equilibrium variables and the transformation matrices are deleted after use. Also, S, sigma, kappa_n and kappa_g can be neglected as they do not contain derivatives and are therefore much smaller. in both routines calc_memory_eq() and calc_memory_x(), however, a 50% safety factor is used to account for this somewhat.| [in] | n_par_x | number of parallel points to be divided |
| [in] | arr_size | array size (using loc_n_r) for eq_2 and X_1 variables |
| [in,out] | n_div | final number of divisions |
| [in] | n_div_max | maximum n_div |
| [in] | n_par_x_base | base n_par_X, undivisible |
| [in] | range_name | name of range |
Definition at line 6564 of file eq_ops.f90.
| integer function, public eq_ops::flux_q_plot | ( | type(grid_type), intent(in) | grid_eq, |
| type(eq_1_type), intent(in) | eq ) |
Plots the flux quantities in the solution grid.
q_saf rot_t pres flux_p flux_t | [in] | grid_eq | normal grid |
| [in] | eq | flux equilibrium variables |
Definition at line 3809 of file eq_ops.f90.
| integer function, public eq_ops::j_plot | ( | type(grid_type), intent(inout) | grid_eq, |
| type(eq_1_type), intent(in) | eq_1, | ||
| type(eq_2_type), intent(in) | eq_2, | ||
| integer, intent(in), optional | rich_lvl, | ||
| logical, intent(in), optional | plot_fluxes, | ||
| real(dp), dimension(:,:,:,:), intent(in), optional | xyz ) |
Plots the current.
If multiple equilibrium parallel jobs, every job does its piece, and the results are joined automatically by plot_HDF5.
The outputs are given in contra- and covariant components and magnitude in multiple coordinate systems, as indicated in calc_vec_comp().
The starting point is the pressure balance
\[ \nabla p = \vec{J} \times \vec{B}, \]
which, using \(\vec{B} = \frac{\vec{e}_{\theta}}{\mathcal{J}}\), reduces to
\[J^\alpha = -p'. \]
Furthermore, the current has to lie in the magnetic flux surfaces:
\[J^\psi = 0. \]
Finally, the parallel current \(\sigma\) gives an expression for the last contravariant component:
\[J^\theta = \frac{\sigma}{\mathcal{J}} + p' \frac{B_\alpha}{B_\theta}. \]
From these, the contravariant components can be calculated as
\[J_i = J^\alpha g_{\alpha,i} + J^\theta g_{\theta,i}. \]
These are then all be transformed to the other coordinate systems.
| [in,out] | grid_eq | equilibrium grid |
| [in] | eq_1 | flux equilibrium variables |
| [in] | eq_2 | metric equilibrium variables |
| [in] | rich_lvl | Richardson level |
| [in] | plot_fluxes | plot the fluxes |
| [in] | xyz | X, Y and Z of grid |
Definition at line 5801 of file eq_ops.f90.
| integer function, public eq_ops::kappa_plot | ( | type(grid_type), intent(inout) | grid_eq, |
| type(eq_1_type), intent(in) | eq_1, | ||
| type(eq_2_type), intent(in) | eq_2, | ||
| integer, intent(in), optional | rich_lvl, | ||
| real(dp), dimension(:,:,:,:), intent(in), optional | xyz ) |
Plots the curvature.
If multiple equilibrium parallel jobs, every job does its piece, and the results are joined automatically by plot_HDF5.
The outputs are given in contra- and covariant components and magnitude in multiple coordinate systems, as indicated in calc_vec_comp().
The starting point is the curvature, given by
\[\vec{\kappa} = \kappa_n \frac{\nabla \psi}{ \left|\nabla \psi\right|^2 } + \kappa_g \frac{\nabla \psi \times \vec{B}}{B^2} , \]
which can be used to find the covariant and contravariant components in Flux coordinates.
These are then transformed to Cartesian coordinates and plotted.
| [in,out] | grid_eq | equilibrium grid |
| [in] | eq_1 | metric equilibrium variables |
| [in] | eq_2 | metric equilibrium variables |
| [in] | rich_lvl | Richardson level |
| [in] | xyz | X, Y and Z of grid |
Definition at line 5969 of file eq_ops.f90.
| subroutine, public eq_ops::normalize_input |
Normalize input quantities.
Definition at line 5642 of file eq_ops.f90.
| logical, public eq_ops::debug_calc_derived_q = .false. |
plot debug information for calc_derived_q()
Definition at line 30 of file eq_ops.f90.
| logical, public eq_ops::debug_create_vmec_input = .false. |
plot debug information for create_vmec_input()
Definition at line 34 of file eq_ops.f90.
| logical, public eq_ops::debug_j_plot = .false. |
plot debug information for j_plot()
Definition at line 32 of file eq_ops.f90.