Numerical utilities related to the grids and different coordinate systems. More...

Interfaces and Types
interface	calc_eqd_grid
	Calculate grid of equidistant points, where optionally the last point can be excluded. More...

interface	calc_tor_diff
	Calculates the toroidal difference for a magnitude calculated on three toroidal points: two extremities and one in the middle. More...

interface	coord_e2f
	Converts Equilibrium coordinates \(\left(r,\theta,\zeta\right)_\text{E}\). to Flux coordinates \(\left(r,\theta,\zeta\right)_\text{F}\). More...

interface	coord_f2e
	Converts Flux coordinates \(\left(r,\theta,\zeta\right)_\text{F}\) to Equilibrium coordinates \(\left(r,\theta,\zeta\right)_\text{E}\). More...

Functions/Subroutines
integer function	coord_e2f_rtz (grid_eq, r_E, theta_E, zeta_E, r_F, theta_F, zeta_F, r_E_array, r_F_array)
	version with r, theta and zeta More...

integer function	coord_e2f_r (grid_eq, r_E, r_F, r_E_array, r_F_array)
	version with only r More...

integer function, public	calc_xyz_grid (grid_eq, grid_XYZ, X, Y, Z, L, R)
	Calculates \(X\), \(Y\) and \(Z\) on a grid `grid_XYZ`, determined through its E(quilibrium) coordinates. More...

integer function, public	extend_grid_f (grid_in, grid_ext, grid_eq, n_theta_plot, n_zeta_plot, lim_theta_plot, lim_zeta_plot)
	Extend a grid angularly. More...

integer function, public	copy_grid (grid_A, grid_B, lims_B, i_lim)
	Copy a grid A to a new grid B, that was not yet initialized. More...

integer function, public	calc_int_vol (ang_1, ang_2, norm, J, f, f_int)
	Calculates volume integral on a 3D grid. More...

integer function, public	trim_grid (grid_in, grid_out, norm_id)
	Trim a grid, removing any overlap between the different regions. More...

integer function, public	untrim_grid (grid_in, grid_out, size_ghost)
	Untrims a trimmed grid by introducing an asymmetric ghost regions at the right. More...

integer function, public	calc_vec_comp (grid, grid_eq, eq_1, eq_2, v_com, norm_disc_prec, v_mag, base_name, max_transf, v_flux_tor, v_flux_pol, XYZ, compare_tor_pos)
	Calculates contra- and covariant components of a vector in multiple coordinate systems. More...

integer function, public	calc_n_par_x_rich (n_par_X_rich, only_half_grid)
	Calculates the local number of parallel grid points for this Richardson level, taking into account that it ould be half the actual number. More...

integer function, public	nufft (x, f, f_F, plot_name)
	calculates the cosine and sine mode numbers of a function defined on a non-regular grid. More...

subroutine, public	find_compr_range (r_F, lim_r, lim_id)
	finds smallest range that comprises a minimum and maximum value. More...

integer function, public	calc_arc_angle (grid, eq_1, eq_2, arc, use_E)
	Calculate arclength angle. More...

Variables
logical, public	debug_calc_int_vol = .false.
	plot debug information for calc_int_vol() More...

logical, public	debug_calc_vec_comp = .false.
	plot debug information for calc_vec_comp() More...

Detailed Description

Numerical utilities related to the grids and different coordinate systems.

Function/Subroutine Documentation

◆ calc_arc_angle()

integer function, public grid_utilities::calc_arc_angle	(	type(grid_type), intent(in)	grid,
		type(eq_1_type), intent(in)	eq_1,
		type(eq_2_type), intent(in)	eq_2,
		real(dp), dimension(:,:,:), intent(inout), allocatable	arc,
		logical, intent(in), optional	use_E
	)

Calculate arclength angle.

This angle is based on the calculation of the arclength from the start of the grid, and then normalizing it from \(\min(\gamma)\ldots\max(\gamma)\), where \(\gamma\) is the parallel angle.

By default, the Flux variables are used, but this can be changed with use_E.

Note: For field-aligned grids the projection is taken on the poloidal cross-section. For non-axisymmetric equilibria, this makes little sense.

Returns: ierr

Parameters

[in]	grid	equilibrium grid
[in]	eq_1	flux equilibrium variables
[in]	eq_2	metric equilibrium variables
[in,out]	arc	arclength angle
[in]	use_e	use E variables instead of F

Definition at line 3039 of file grid_utilities.f90.

Here is the call graph for this function:

◆ calc_int_vol()

integer function, public grid_utilities::calc_int_vol	(	real(dp), dimension(:,:,:), intent(in)	ang_1,
		real(dp), dimension(:,:,:), intent(in)	ang_2,
		real(dp), dimension(:), intent(in)	norm,
		real(dp), dimension(:,:,:), intent(in)	J,
		complex(dp), dimension(:,:,:,:), intent(in)	f,
		complex(dp), dimension(:), intent(inout)	f_int
	)

Calculates volume integral on a 3D grid.

Two angular and one normal variable have to be provided on a the grid.

If the i'th dimension of the grid is equal to one, the function to be integrated is assumed not to vary in this dimension.

Furthermore, if i is 1 or 2, the corresponding i'th (angular) variable is the only variable that is assumed to vary in that dimension. The other angular variable as well as the normal variable are assumed to be constant like the function itself, resulting in a factor \(2 \pi\). However, if i is 3, an error is displayed as this does not represent a physical situation.

A common case through which to understand this is the axisymmetric case where the first angular variable \(\theta\) varies in the dimensions 1 and 3, the second angular variable \(\zeta\) varies only in the dimension 2, and the normal variable only varies in the dimension 3.

Alternatively, there is the case of a grid-aligned set of coordinates \(\theta\) and \(\alpha\), where the first dimension corresponds to the direction along the magnetic field line, the second to the geodesical direction and the third to the normal direction. If the calculations for different field lines are decoupled, the variation in the second dimension is not taken into account and no integration happens along it.

Internally, the angular variables and the normal variable are related to the coordinates \(\left(x,y,z\right)\) that correspond to the three dimensions. They thus form a computational orthogonal grid to which the original coordinates are related through the transformation of Jacobians:

\[\mathcal{J}_{xyz} = \mathcal{J}_\text{F} \frac{\partial r_\text{F}}{\partial z} \left(\frac{\partial \gamma_1}{\partial x} \frac{\gamma_2}{\partial y}- \frac{\partial \gamma_1}{\partial y} \frac{\gamma_2}{\partial x}\right) , \]

where \(\gamma_1\) and \(\gamma_2\) are ang_1 and ang_2, so that the integral becomes

\[ \sum_{x,y,z} \left[ f\left(x,y,z\right) \mathcal{J}_\text{F} \frac{\partial r_\text{F}}{\partial z} \left(\frac{\partial \gamma_1}{\partial x} \frac{\gamma_2}{\partial y}- \frac{\partial \gamma_1}{\partial y} \frac{\gamma_2}{\partial x}\right) \Delta_x \Delta_y \Delta_z \right] \]

where \(\Delta_x\), \(\Delta_y\) and \(\Delta_z\) are all trivially equal to 1.

The integrand has to be evaluated at the intermediate positions inside the cells. This is done by taking the average of the \(2^3=8\) points for fj = \(f J_\text{F}\) as well as the transformation of the Jacobian to \(\left(x,y,z\right)\) coordinates.

See also: See grid_vars.grid_type for a discussion on ang_1 and ang_2.

Note

if the coordinates are independent, this method is equivalent to the repeated numerical integration using the trapezoidal method, NOT Simpson's 3/8 rule!
by setting debug_calc_int_vol, this method can be compared to the trapezoidal and simple method for independent coordinates, again NOT for Simpson's 3/8 rule!
The Simpson's 3/8 rule could be developed but it is not of great importance.

Returns: ierr

Parameters

[in]	ang_1	coordinate variable \(\gamma_1\)
[in]	ang_2	coordinate variable \(\gamma_2\)
[in]	norm	coordinate variable \(r_\text{F}\)
[in]	j	Jacobian
[in]	f	input f(n_par,n_geo,n_r,size_X^2)
[in,out]	f_int	output integrated f

Definition at line 1320 of file grid_utilities.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ calc_n_par_x_rich()

integer function, public grid_utilities::calc_n_par_x_rich	(	integer, intent(inout)	n_par_X_rich,
		logical, intent(in), optional	only_half_grid
	)

Calculates the local number of parallel grid points for this Richardson level, taking into account that it ould be half the actual number.

See also: calc_ang_grid_eq_b().

Returns: ierr

Parameters

[in,out]	n_par_x_rich	n_par_X for this Richardson level
[in]	only_half_grid	calculate only half grid with even points

Definition at line 2862 of file grid_utilities.f90.

Here is the caller graph for this function:

◆ calc_vec_comp()

integer function, public grid_utilities::calc_vec_comp	(	type(grid_type), intent(in)	grid,
		type(grid_type), intent(in)	grid_eq,
		type(eq_1_type), intent(in)	eq_1,
		type(eq_2_type), intent(in)	eq_2,
		real(dp), dimension(:,:,:,:,:), intent(inout)	v_com,
		integer, intent(in)	norm_disc_prec,
		real(dp), dimension(:,:,:), intent(inout), optional	v_mag,
		character(len=*), intent(in), optional	base_name,
		integer, intent(in), optional	max_transf,
		real(dp), dimension(:,:), intent(inout), optional, allocatable	v_flux_tor,
		real(dp), dimension(:,:), intent(inout), optional, allocatable	v_flux_pol,
		real(dp), dimension(:,:,:,:), intent(in), optional	XYZ,
		logical, intent(in), optional	compare_tor_pos
	)

Calculates contra- and covariant components of a vector in multiple coordinate systems.

Chain of coordinate systems considered:

Flux F: \((\alpha,\psi,\theta)\),
Magnetic M: \((\psi,\theta,\zeta)\),
- H: \((\psi,\theta,\phi)\) if HELENA is used,
- V: \((r,\theta,\zeta)\) if VMEC is used,
Cylindrical C: \((R,\varphi,Z)\),
Cartesian X: \((X,Y,Z)\), as well as the magnitude, starting from the input in Flux coordinates. Note that the Cartesian components in X, Y and Z can be used to plot the real vector and that covariant components should be equal to contravariant components in this coordinate system.

Also, the fluxes can be calculated and plot by providing a base_name.

By default, the cartesian components are returned, but this can be indicated differently by providing max_transf.

Note

Plots for different Richardson levels can be combined to show the total grid by just plotting them all individually and sequentially.
The metric factors and transformation matrices have to be allocated. They can be calculated using the routines from eq_ops, for deriv = [0,0,0].
For VMEC, the trigonometric factors of grid_XYZ must be calculated beforehand. Optionally, by providing X, Y and Z, the ones calculated in this routine are overwritten. This is useful for, for example, slab geometries.
The normalization factors are taken into account and the output is transformed back to unnormalized values.

Returns: ierr

Parameters

[in]	grid	grid on which vector is calculated
[in]	grid_eq	grid on which equilibrium variables are calculated
[in]	eq_1	flux equilibrium quantities
[in]	eq_2	metric equilibrium variables
[in,out]	v_com	covariant and contravariant components `(dim1,dim2,dim3,3,2)`
[in]	norm_disc_prec	precision for normal derivatives
[in,out]	v_mag	magnitude `(dim1,dim2,dim3)`
[in]	base_name	base name for output plotting
[in]	max_transf	maximum transformation level (2: Magnetic, 3: Equilibrium, 4: Cylindrical, 5: Cartesian [def])
[in,out]	v_flux_pol	poloidal flux as function of normal coordinate for all toroidal positions
[in,out]	v_flux_tor	toroidal flux as function of normal coordinate for all poloidal positions
[in]	xyz	\(X\), \(Y\) and \(Z\) of grid
[in]	compare_tor_pos	compare toroidal positions

Definition at line 1859 of file grid_utilities.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ calc_xyz_grid()

integer function, public grid_utilities::calc_xyz_grid	(	type(grid_type), intent(in)	grid_eq,
		type(grid_type), intent(in)	grid_XYZ,
		real(dp), dimension(:,:,:), intent(inout)	X,
		real(dp), dimension(:,:,:), intent(inout)	Y,
		real(dp), dimension(:,:,:), intent(inout)	Z,
		real(dp), dimension(:,:,:), intent(inout), optional	L,
		real(dp), dimension(:,:,:), intent(inout), optional	R
	)

Calculates \(X\), \(Y\) and \(Z\) on a grid grid_XYZ, determined through its E(quilibrium) coordinates.

Furthermore, a grid grid_eq must be provided, which is the grid in which the variables concerning \(R\) and \(Z\) are tabulated, i.e. the full equilibrium grid in E(quilibrium) coordinates.

Of this grid, however, only r_E is used, and the rest ignored. It can therefore be provided without the angular part, i.e. by reconstructing it with a subset.

If VMEC is the equilibrium model, this routine also optionally calculates \(\lambda\) on the grid, as this is also needed some times. for HELENA this variable is not used.

Optionally, \(R\) and \(\lambda\) (only for VMEC) can be returned.

Note

For VMEC, the trigonometric factors of grid_XYZ must be calculated beforehand.
The normalization factor R_0 for length is taken into account and the output is transformed back to unnormalized values.

Returns: ierr

Parameters

[in]	grid_eq	equilibrium grid
[in]	grid_xyz	grid for which to calculate \(X\), \(Y\), \(Z\) and optionally \(L\)
[in,out]	x	\(X\) of grid
[in,out]	y	\(Y\) of grid
[in,out]	z	\(Z\) of grid
[in,out]	l	\(\lambda\) of grid
[in,out]	r	\(R\) of grid

Definition at line 799 of file grid_utilities.f90.

Here is the caller graph for this function:

◆ coord_e2f_r()

integer function grid_utilities::coord_e2f_r	(	type(grid_type), intent(in)	grid_eq,
		real(dp), dimension(:), intent(in)	r_E,
		real(dp), dimension(:), intent(inout)	r_F,
		real(dp), dimension(:), intent(in), optional, target	r_E_array,
		real(dp), dimension(:), intent(in), optional, target	r_F_array
	)

version with only r

Parameters

[in]	grid_eq	equilibrium grid (for normal local limits)
[in]	r_e	(local) \(r_\text{E}\)
[in,out]	r_f	(local) \(r_\text{F}\)
[in]	r_e_array	optional array that defines mapping between two coord. systems
[in]	r_f_array	optional array that defines mapping between two coord. systems

Definition at line 473 of file grid_utilities.f90.

Here is the caller graph for this function:

◆ coord_e2f_rtz()

integer function grid_utilities::coord_e2f_rtz	(	type(grid_type), intent(in)	grid_eq,
		real(dp), dimension(:), intent(in)	r_E,
		real(dp), dimension(:,:,:), intent(in)	theta_E,
		real(dp), dimension(:,:,:), intent(in)	zeta_E,
		real(dp), dimension(:), intent(inout)	r_F,
		real(dp), dimension(:,:,:), intent(inout)	theta_F,
		real(dp), dimension(:,:,:), intent(inout)	zeta_F,
		real(dp), dimension(:), intent(in), optional, target	r_E_array,
		real(dp), dimension(:), intent(in), optional, target	r_F_array
	)

version with r, theta and zeta

Parameters

[in]	grid_eq	equilibrium grid (for normal local limits)
[in]	r_e	(local) \(r_\text{E}\)
[in]	theta_e	\(\theta_\text{E}\)
[in]	zeta_e	\(\zeta_\text{E}\)
[in,out]	r_f	(local) \(r_\text{F}\)
[in,out]	theta_f	\(\theta_\text{F}\)
[in,out]	zeta_f	\(\zeta_\text{F}\)
[in]	r_e_array	optional array that defines mapping between two coord. systems
[in]	r_f_array	optional array that defines mapping between two coord. systems

Definition at line 359 of file grid_utilities.f90.

Here is the call graph for this function:

◆ copy_grid()

integer function, public grid_utilities::copy_grid	(	class(grid_type), intent(in)	grid_A,
		class(grid_type), intent(inout)	grid_B,
		integer, dimension(3,2), intent(in), optional	lims_B,
		integer, dimension(2), intent(in), optional	i_lim
	)

Copy a grid A to a new grid B, that was not yet initialized.

This new grid can contain a subsection of the previous grid in all dimensions. It can also be a divided grid, by providing the limits

Note

The normal limits for the divided grid should be given in terms of the normal dimension of grid B.
if the grids are not purely normal, the procedure can currently only handle the copying of a grid_A that is not divided.

Returns: ierr

Parameters

[in]	grid_a	grid to be initialized
[in,out]	grid_b	grid to be initialized
[in]	lims_b	ranges for subset of grid
[in]	i_lim	min. and max. local normal index

Definition at line 1189 of file grid_utilities.f90.

Here is the caller graph for this function:

◆ extend_grid_f()

integer function, public grid_utilities::extend_grid_f	(	type(grid_type), intent(in)	grid_in,
		type(grid_type), intent(inout)	grid_ext,
		type(grid_type), intent(in), optional	grid_eq,
		integer, intent(in), optional	n_theta_plot,
		integer, intent(in), optional	n_zeta_plot,
		real(dp), dimension(2), intent(in), optional	lim_theta_plot,
		real(dp), dimension(2), intent(in), optional	lim_zeta_plot
	)

Extend a grid angularly.

This is done using equidistant variables of n_theta_plot and n_zeta_plot angular and own loc_n_r points in F coordinates.

Optionally, the grid can also be converted to E coordinates if the equilibrium grid is provided. This is required for many operations, such as the calculation of \(X\), \(Y\) and \(Z\) through calc_XYZ_grid().

Note: For VMEC, it can be slow, as grid_utilities.coord_f2e() is used.

Returns: ierr

Parameters

[in]	grid_in	grid to be extended
[in,out]	grid_ext	extended grid
[in]	grid_eq	equilibrium grid
[in]	n_theta_plot	number of poins in theta direction
[in]	n_zeta_plot	number of poins in zeta direction
[in]	lim_theta_plot	limits in theta
[in]	lim_zeta_plot	limits in zeta

Definition at line 1096 of file grid_utilities.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ find_compr_range()

subroutine, public grid_utilities::find_compr_range	(	real(dp), dimension(:), intent(in)	r_F,
		real(dp), dimension(2), intent(in)	lim_r,
		integer, dimension(2), intent(inout)	lim_id
	)

finds smallest range that comprises a minimum and maximum value.

Parameters

[in]	r_f	all values of coordinate
[in]	lim_r	limiting range
[in,out]	lim_id	limiting indices

Definition at line 2995 of file grid_utilities.f90.

Here is the caller graph for this function:

◆ nufft()

integer function, public grid_utilities::nufft	(	real(dp), dimension(:), intent(in)	x,
		real(dp), dimension(:), intent(in)	f,
		real(dp), dimension(:,:), intent(inout), allocatable	f_F,
		character(len=*), intent(in), optional	plot_name
	)

calculates the cosine and sine mode numbers of a function defined on a non-regular grid.

If a plot name is provided, the modes are plotted.

Note: The fundamental interval is assumed to be \(0\ldots 2\pi\) but there should be no overlap between the first and last point.

Returns: ierr

Parameters

[in]	x	coordinate values
[in]	f	function values
[in,out]	f_f	Fourier modes for \(\cos\) and \(\sin\)
[in]	plot_name	name of possible plot

Definition at line 2901 of file grid_utilities.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ trim_grid()

integer function, public grid_utilities::trim_grid	(	type(grid_type), intent(in)	grid_in,
		type(grid_type), intent(inout)	grid_out,
		integer, dimension(2), intent(inout), optional	norm_id
	)

Trim a grid, removing any overlap between the different regions.

by default, the routine assumes a symmetric ghost region and cuts as many grid points from the end of the previous process as from the beginning of the next process, but if the number of overlapping grid points is odd, the next process looses one more point.

optionally, the trimmed indices in the normal direction can be provided in norm_id, i.e. the indices in the old, untrimmed grid that correspond to the start and end indices of the trimmed grid. E.g. if

proc 0: 3 ... 25
proc 1: 20 ... 50

then the trimmed grid will be:

proc 0: 3 ... 22
proc 1: 23 ... 50

which is shifted down by 2 to

proc 0: 1 ... 20
proc 1: 21 ... 48

in the trimmed grid. The indices of the previous step (3 & 22 and 23 & 50) are saved in norm_id.

Returns: ierr

Parameters

[in]	grid_in	input grid
[in,out]	grid_out	trimmed grid
[in,out]	norm_id	normal indices corresponding to trimmed part

Definition at line 1636 of file grid_utilities.f90.

Here is the caller graph for this function:

◆ untrim_grid()

integer function, public grid_utilities::untrim_grid	(	type(grid_type), intent(in)	grid_in,
		type(grid_type), intent(inout)	grid_out,
		integer, intent(in)	size_ghost
	)

Untrims a trimmed grid by introducing an asymmetric ghost regions at the right.

The width of the ghost region has to be provided.

Note

the ghosted grid should be deallocated (with dealloc_grid()).
the input grid has to be trimmed.

Returns: ierr

Parameters

[in]	grid_in	input grid
[in,out]	grid_out	ghosted grid
[in]	size_ghost	width of ghost region

Definition at line 1764 of file grid_utilities.f90.

Here is the caller graph for this function:

Variable Documentation

◆ debug_calc_int_vol

logical, public grid_utilities::debug_calc_int_vol = .false.

plot debug information for calc_int_vol()

Note: Debug version only

Definition at line 25 of file grid_utilities.f90.

◆ debug_calc_vec_comp

logical, public grid_utilities::debug_calc_vec_comp = .false.

plot debug information for calc_vec_comp()

Note: Debug version only

Definition at line 27 of file grid_utilities.f90.

Interfaces and Types

Functions/Subroutines

Variables

Detailed Description

Function/Subroutine Documentation

◆ calc_arc_angle()

◆ calc_int_vol()

◆ calc_n_par_x_rich()

◆ calc_vec_comp()

◆ calc_xyz_grid()

◆ coord_e2f_r()

◆ coord_e2f_rtz()

◆ copy_grid()

◆ extend_grid_f()

◆ find_compr_range()

◆ nufft()

◆ trim_grid()

◆ untrim_grid()

Variable Documentation

◆ debug_calc_int_vol

◆ debug_calc_vec_comp