eq_ops Module Reference

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	calc_g_c
	Calculate the lower metric elements in the C(ylindrical) coordinate system. More...
interface	calc_g_f
	Calculate the metric coefficients in the F(lux) coordinate system. More...
interface	calc_g_h
	Calculate the lower metric coefficients in the equilibrium H(ELENA) coordinate system. More...
interface	calc_g_v
	Calculate the lower metric coefficients in the equilibrium V(MEC) coordinate system. More...
interface	calc_jac_f
	Calculate \(\mathcal{J}_\text{F}\), the jacobian in Flux coordinates. More...
interface	calc_jac_h
	Calculate \(\mathcal{J}_\text{H}\), the jacobian in HELENA coordinates. More...
interface	calc_jac_v
	Calculate \(\mathcal{J}_\text{V}\), the jacobian in V(MEC) coordinates. More...
interface	calc_rzl
	Calculate \(R\), \(Z\) & \(\lambda\) and derivatives in VMEC coordinates. 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...
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	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...

Functions/Subroutines
integer function	create_vmec_input (grid_eq, eq_1)
	Creates a VMEC input file. More...
integer function, public	flux_q_plot (grid_eq, eq)
	Plots the flux quantities in the solution grid. More...
integer function, public	calc_derived_q (grid_eq, eq_1, eq_2)
	Calculates derived equilibrium quantities system. More...
integer function, public	calc_normalization_const ()
	Sets up normalization constants. More...
subroutine, public	normalize_input ()
	Normalize input quantities. More...
integer function, public	b_plot (grid_eq, eq_1, eq_2, rich_lvl, plot_fluxes, XYZ)
	Plots the magnetic fields. More...
integer function, public	j_plot (grid_eq, eq_1, eq_2, rich_lvl, plot_fluxes, XYZ)
	Plots the current. More...
integer function, public	kappa_plot (grid_eq, eq_1, eq_2, rich_lvl, XYZ)
	Plots the curvature. More...
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. More...
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. More...
integer function, public	calc_eq_jobs_lims (n_par_X, n_div)
	Calculate eq_jobs_lims. More...
integer function	test_t_ef (grid_eq, eq_1, eq_2)
	See if T_EF it complies with the theory of [17]. More...
integer function	test_d12h_h (grid_eq, eq)
	Tests whether \( \frac{\partial^2}{\partial u_i \partial u_j} h_\text{H} \) is calculated correctly. More...
integer function	test_jac_f (grid_eq, eq_1, eq_2)
	Performs tests on \( \mathcal{J}_\text{F}\). More...
integer function	test_g_v (grid_eq, eq)
	Tests whether \(g_\text{V}\) is calculated correctly. More...
integer function	test_jac_v (grid_eq, eq)
	Tests whether \(\mathcal{J}_\text{V}\) is calculated correctly. More...
integer function	test_b_f (grid_eq, eq_1, eq_2)
	Tests whether \(\vec{B}_\text{F}\) is calculated correctly. More...
integer function	test_p (grid_eq, eq_1, eq_2)
	Performs tests on pressure balance. More...

Variables
logical, public	debug_calc_derived_q = .false.
	plot debug information for calc_derived_q() More...
logical, public	debug_j_plot = .false.
	plot debug information for j_plot() More...
logical, public	debug_create_vmec_input = .false.
	plot debug information for create_vmec_input() More...

Detailed Description

Operations on the equilibrium variables.

Function/Subroutine Documentation

b_plot()

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.

Note

1. Vector plots for different Richardson levels can be combined to show the total grid by just plotting them all individually.
2. The metric factors and transformation matrices have to be allocated.

Returns

ierr

Parameters

[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 5700 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

calc_derived_q()

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.

• magnetic shear S
• normal curvature kappa_n
• geodesic curvature kappa_g
• parallel current 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.

VMEC

The shear

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.

The curvature

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

1. Set up \(\mathcal{D}\vec{e}_\text{E}\), i.e. as a function of the three cylindrical covariant basis vectors.
2. Set up correctcion by the derivatives of the transformation matrix, with a single summation.
3. Decompose the normal \(\frac{\nabla \psi}{\left|\nabla\psi\right|^2}\) and geodesic vectors \(\frac{\nabla \psi \times \vec{B}}{B^2}\) as a function of the cylindrical contravariant basis vectors.
4. Double summation to reduce term ~ \(\mathcal{D}\vec{e}_\text{E}\) and correction by single summation to reduce term ~ \(\mathcal{D}\mathcal{T}\).
5. Dot these for each of the four (actually three because of symmetry) elements of \(\mathcal{D}\vec{e}_\text{E}\).
6. Divide by \(\left|\vec{e}_\theta\right|^2\).
7. Possibly correct for toroidal flux.

An advantage of using this formulation is that no normal derivatives are needed, so that nothing has to implicitely cancel out.

The parallel current

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.

HELENA

The parallel current

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 shear

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.

The curvature

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} \)

Note

1. If the toroidal flux is used instead, the actual curvature obviously is unchanged, which implies that the normal component has to be multiplied by the safety factor and the geodesic component has to be divided by it.
2. The formulas for the normal and geodesic basis vectors for VMEC are
   • \(\frac{\nabla \psi}{\left|\nabla\psi\right|^2} = -\frac{\mathcal{J} \left(1 + \lambda_\theta\right)} {Z_\theta^2 + R_\theta^2 + \left(R_\zeta Z_\theta - R_\theta Z_\zeta \right)/R^2} \frac{1}{R} \left( -Z_\theta , R_\zeta Z_\theta - R_\theta Z_\zeta, R_\theta\right)_\text{C}\)
   • \(\frac{\nabla \psi \times \vec{B}}{B^2} = \frac{1}{g_{\theta\theta}} \left(P R_\theta - Q R_\zeta, -Q R^2, P Z_\theta - Q Z_\zeta\right)_\text{C}\)
   • \(Q = g_{\theta\theta} - q g_{\theta\alpha}\) (using right-handed flux \(q_\text{F}\), which is the inverse of the left-handed VMEC \(q_\text{V}\))
   • \(P = \frac{Q \lambda_\zeta - g_{\alpha\theta}}{1+\lambda_\theta}\) (where the derivatives of \(R\), \(Z\) and \(\lambda\) are in VMEC coordinates).
3. The formulas for the normal and geodesic basis vectors for HELENA are
   • \(\frac{\nabla \psi}{\left|\nabla\psi\right|^2} = -\frac{1}{Z_\theta^2 + R_\theta^2} \frac{q R}{F} \left( -Z_\theta , 0, R_\theta\right)_\text{C}\)
   • \(\frac{\nabla \psi \times \vec{B}}{B^2} = \frac{q R^2}{R_\theta^2 + Z_\theta^2 + q^2 R^2} \left(-R_\theta, -\frac{R_\theta^2 + Z_\theta^2}{q}, -Z_\theta\right)_\text{C}\)

Parameters

[in]	grid_eq	equilibrium grid
[in]	eq_1	flux equilibrium variables
[in,out]	eq_2	metric equilibrium variables

Definition at line 4287 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

calc_eq_jobs_lims()

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():

• for magn_int_style=1 (trapezoidal), this is 1 point,
• for 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.

Note

The 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.

Returns

ierr

Parameters

[in]	n_par_x	number of parallel points in this Richardson level
[in]	n_div	nr. of divisions of parallel ranges

Definition at line 6701 of file eq_ops.f90.

Here is the caller graph for this function:

calc_normalization_const()

integer function, public eq_ops::calc_normalization_const

Sets up normalization constants.

• VMEC version (eq_style=1):
  Normalization depends on norm_style:
  1. MISHKA
     • R_0: major radius (= average magnetic axis)
     • B_0: B on magnetic axis (theta = zeta = 0)
     • pres_0: reference pressure (= B_0^2/mu_0)
     • psi_0: reference flux (= R_0^2 B_0)
     • rho_0: reference mass density
  2. COBRA
     • R_0: major radius (= average geometric axis)
     • pres_0: pressure on magnetic axis
     • B_0: reference magnetic field (= sqrt(2pres_0mu_0 / beta))
     • psi_0: reference flux (= R_0^2 B_0 / aspr^2)
     • rho_0: reference mass density where aspr (aspect ratio) and beta are given by VMEC.
• HELENA version (eq_style=2):
  MISHKA Normalization is used by default and does not depend on norm_style

See also

read_hel()

Note

rho_0 is not given through by the equilibrium codes and should be user-supplied

Returns

ierr

Definition at line 5405 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

create_vmec_input()

integer function eq_ops::create_vmec_input	(	type(grid_type), intent(in)	grid_eq,
		type(eq_1_type), intent(in)	eq_1
	)

Creates a VMEC input file.

Optionally, a perturbation can be added: Either the displacement of the plasma position can be described (pert_style 1), or ripple in the toroidal magnetic field (pert_style 2), with a fixed toroidal mode number.

Both perturbation styles can have various prescription types:

1. file with Fourier modes in the geometrical angular coordinate
2. same but manually
3. file with perturbation from a 2-D map in the geometric angular coordinate.

For pert_style 2, a file has to be provided that describes the translation between position perturbation and magnetic perturbation for curves of constant geometrical angle. This file can be generated for an already existing ripple case using POST with –compare_tor_pos with n_zeta_plot = 3 and min_theta_plot and max_theta_plot indicating half a ripple period.

The output from this VMEC run can then be used to iteratively create a new file to translate toroidal magnetic field ripple to position perturbation.

Note

Meaning of the indices of B_F, B_F_dum:

• (pol modes, cos/sin) for B_F_dum
• (tor_modes, pol modes, cos/sin (m theta), R/Z) for B_F

Parameters

[in]	grid_eq	equilibrium grid varibles
[in]	eq_1	flux equilibrium quantities

Definition at line 784 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

delta_r_plot()

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.

Note

The metric factors and transformation matrices have to be allocated.

Returns

ierr

Parameters

[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 6173 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

divide_eq_jobs()

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:

• for magn_int_style=1 (trapezoidal), this is 1 point,
• for 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:

• If the first job of the parallel jobs and not the first Richardson level: add half of the integrated tensorial perturbation quantities of the previous level.
• If not the first job of the parallel jobs, add the integrated tensorial perturbation quantities to those of the previous parallel job, same Richardson level.

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.

Note

For PB3D, only the variables 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.

Returns

ierr

Parameters

[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 6566 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

flux_q_plot()

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.

• safety factor q_saf
• rotational transform rot_t
• pressure pres
• poloidal flux flux_p
• toroidal flux flux_t

Returns

ierr

Parameters

[in]	grid_eq	normal grid
[in]	eq	flux equilibrium variables

Definition at line 3810 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

j_plot()

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.

Note

1. Vector plots for different Richardson levels can be combined to show the total grid by just plotting them all individually.
2. The metric factors and transformation matrices have to be allocated.

Returns

ierr

Parameters

[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 5803 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

kappa_plot()

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.

Note

1. Vector plots for different Richardson levels can be combined to show the total grid by just plotting them all individually.
2. The metric factors and transformation matrices have to be allocated.

Returns

ierr

Parameters

[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 5971 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

normalize_input()

subroutine, public eq_ops::normalize_input

Normalize input quantities.

See also

calc_normalization_const()

Definition at line 5643 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

test_b_f()

integer function eq_ops::test_b_f	(	type(grid_type), intent(in)	grid_eq,
		type(eq_1_type), intent(in)	eq_1,
		type(eq_2_type), intent(in)	eq_2
	)

Tests whether \(\vec{B}_\text{F}\) is calculated correctly.

Note

Debug version only

Returns

ierr

Parameters

[in]	grid_eq	equilibrium grid
[in]	eq_1	flux equilibrium
[in]	eq_2	metric equilibrium

Definition at line 7344 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

test_d12h_h()

integer function eq_ops::test_d12h_h	(	type(grid_type), intent(in)	grid_eq,
		type(eq_2_type), intent(in)	eq
	)

Tests whether \( \frac{\partial^2}{\partial u_i \partial u_j} h_\text{H} \) is calculated correctly.

Note

Debug version only

Returns

ierr

Parameters

[in]	grid_eq	equilibrium grid
[in]	eq	metric equilibrium

Definition at line 6953 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

test_g_v()

integer function eq_ops::test_g_v	(	type(grid_type), intent(in)	grid_eq,
		type(eq_2_type), intent(in)	eq
	)

Tests whether \(g_\text{V}\) is calculated correctly.

Note

Debug version only

Returns

ierr

Parameters

[in]	grid_eq	equilibrium grid
[in]	eq	metric equilibrium

Definition at line 7181 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

test_jac_f()

integer function eq_ops::test_jac_f	(	type(grid_type), intent(in)	grid_eq,
		type(eq_1_type), intent(in), target	eq_1,
		type(eq_2_type), intent(in)	eq_2
	)

Performs tests on \( \mathcal{J}_\text{F}\).

• comparing it with the determinant of \(g_\text{F}\)
• comparing it with the direct formula

Note

Debug version only

Returns

ierr

Parameters

[in]	grid_eq	equilibrium grid
[in]	eq_1	flux equilibrium
[in]	eq_2	metric equilibrium

Definition at line 7064 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

test_jac_v()

integer function eq_ops::test_jac_v	(	type(grid_type), intent(in)	grid_eq,
		type(eq_2_type), intent(in)	eq
	)

Tests whether \(\mathcal{J}_\text{V}\) is calculated correctly.

Note

Debug version only

Returns

ierr

Parameters

[in]	grid_eq	equilibrium grid
[in]	eq	metric equilibrium

Definition at line 7280 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

test_p()

integer function eq_ops::test_p	(	type(grid_type), intent(in)	grid_eq,
		type(eq_1_type), intent(in)	eq_1,
		type(eq_2_type), intent(in)	eq_2
	)

Performs tests on pressure balance.

\[\mu_0 \frac{\partial p}{\partial u^2} = \frac{1}{\mathcal{J}} \left(\frac{\partial B_2}{\partial u^3} - \frac{\partial B_3}{\partial u^2}\right)\]

\[\mu_0 \mathcal{J} \frac{\partial p}{\partial u^3} = 0 \rightarrow \left(\frac{\partial B_1}{\partial u^3} = \frac{\partial B_3}{\partial u^1}\right), \]

working in the (modified) Flux coordinates \(\left(\alpha,\psi,\theta\right)_\text{F}\)

Note

Debug version only

Returns

ierr

Parameters

[in]	grid_eq	equilibrium grid
[in]	eq_1	flux equilibrium variables
[in]	eq_2	metric equilibrium variables

Definition at line 7502 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

test_t_ef()

integer function eq_ops::test_t_ef	(	type(grid_type), intent(in)	grid_eq,
		type(eq_1_type), intent(in)	eq_1,
		type(eq_2_type), intent(in)	eq_2
	)

See if T_EF it complies with the theory of [17].

Note

Debug version only

Returns

ierr

Parameters

[in]	grid_eq	equilibrium grid
[in]	eq_1	flux equilibrium
[in]	eq_2	metric equilibrium

Definition at line 6779 of file eq_ops.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

debug_calc_derived_q

logical, public eq_ops::debug_calc_derived_q = .false.

plot debug information for calc_derived_q()

Note

Debug version only

Definition at line 30 of file eq_ops.f90.

debug_create_vmec_input

logical, public eq_ops::debug_create_vmec_input = .false.

plot debug information for create_vmec_input()

Note

Debug version only

Definition at line 34 of file eq_ops.f90.

debug_j_plot

logical, public eq_ops::debug_j_plot = .false.

plot debug information for j_plot()

Note

Debug version only

Definition at line 32 of file eq_ops.f90.