Procedures

Terminate program correctly if an error has occurred (important in MPI mode!). Uses a MPI_ABORT which terminates FLUXO if a single proc calls this routine.

Add boundary perturbation

Determines if two REAL(wp) numbers are equal up to a specified tolerance (=PP_RealTolerance, normaly set to machine precision) Takes into account that x,y are located in-between [-1;1] Based on Algorithm 139, Kopriva

Applies strong boundary condition to force DOF

Apply preconditioner matrix for single mode of one variable

Computes barycentric (interpolation) weights for interpolation polynomial given by set of nodes. (Algorithm 30, Kopriva book)

change basis from old input base to new base,

compare self and input type base

copy from input type base to self

evaluate all degrees of freedom at all Gauss Points (deriv=0 solution, deriv=1 first derivative d/ds)

evaluate all degrees of freedom at all Gauss Points , for deriv=0 in s and deriv=0,theta,zeta in f

evaluate all degrees of freedom at given s theta zeta position (deriv=0 solution, deriv=1 first derivative d/ds)

finalize the type base

allocate and initialize the type base

test base variable

convert from interpolation points X=> X1_b, Y=> X2_b to fourier modes, given from the input fbase convert to maximum allowable number of modes (ntheta>=2m_max+1, nzeta>=2n_max+1) the final m_max/n_max can be smaller or larger. If larger, a change of base is necessary

initialize class: read file and save data to class structure

Allocate class and call init

get the exact spline antiderivative, with respect to rho2 the knotspan is increased by an extra multiplicity on both ends, and the new coefficients are computed as beta(i) = beta(i-1) + alpha(i)*(t(i+degree+1)-t(i))/(degree+1) From deBoor, "A practical guide to Splines", p.128

evaluate the n-th derivative of the bsplProfile at position s

finalize the type rProfile

initialize the rProfile of type bspline

Build a 1D Vandermonde matrix from an orthonormal Legendre basis to a nodal basis and reverse

Build preconditioner matrices for X1,X2,LA and factorize, for all modes the matrix is only radially dependent, and has the form K_ij = int(s,0,1) d/ds sbase_i(s) (s) d/ds sbase_j(s) + sbase_i(s) ((s) + |Phi'(s)|^2 (-m^2 (s) - n^2 (s) ) ) sbase_j(s) where < > denote an average over the angular coordinates

Builds X1 and X2 in SFL coordinates

Compute Chebychev-Gauss-Lobatto nodes and integration weights (algorithm 27, Kopriva book)

Compute Chebychev-Gauss nodes and integration weights (algorithm 27, Kopriva book)

Check that the TNB frame really has the field periodicity of NFP: assumption for now is that the origin is fixed at rot_origin=(/0.,0.,0./) and the rotation axis is fixed at rot_axis=(/0.,0.,1./) sign of the rotation 'sgn_rot' is now accounted for in the transformation to xhat, so it has to be passed here.

Sample axis and check for zero (<1.e-12) curvature

Compute Clenshaw-Curtis nodes and integration weights

Counts all occurances of string named "key" from inifile and store in "CNTSTR". If keyword "Key" is not found in ini file, the default value "Proposal" is used for "CNTSTR" (error if "Proposal" not given). Inifile was read in before and is stored as list of character strings starting with "FirstString".

Get logical, integer, real, integer array or real array and transform it to string in the proposal format

computes the cross product of to 3 dimensional vectors: cross=v1 x v2

evaluate the n-th derivative of the bsplProfile at position s

finalize the type rProfile

Interpolation of function values f(x_i)=f_i, i=1,n with a cubic spline, given left and right boundary condition types of boundary conditions: 0: not-a-knot 1: f'(x_boundary)=0 2: f''(x_boundary)=0

Remove string "Str" from list of strings witFirstString,h first element "DirstString" and delete string.

compute determinant of 3x3 matrix

Determines if two REAL(wp) numbers are equal up to a given tolerance. Routine requires: x,y > tolerance

@brief Write error/warning message to a given unit. @details

evalute monomial polynomial c_1+c_2x+c_3x^2 ...

evalute first derivative monomial polynomial (c_1+c_2x+c_3x^2) -> (c_2+2c_3x+3c_4x^2 ...

evaluate 1d fourier series from given cos/sin coefficients and mode numbers xn SUM(xc(0:n_max)COS(xn(0:n_max)zeta)+xs(0:n_max)SIN(xn(0:n_max)zeta) evaluate all derivatives 1,2,3 alongside

Evaluate iota profile from restart data

Evaluate phi profile from restart data

Evaluate phiPrime profile from restart data

Evaluate pressure profile from restart data

Evaluate profile/profile derivative given by index "iprof" from restart data

Evaluate auxiliary variables at input state, writes onto module variables!!!

Evaluate 3D MHD energy NOTE: set callEvalaux >0 if not called before for the same Uin !!

Evaluate the variation of the Energy = Force... F_j=-(D_U W(U))_j = -DW(u_h)*testfunc_j NOTE: set callEvalaux TRUE if not called before for the same Uin !!

Evaluate total volume and average surface

Evaluate the basis and all derivatives for a list of (rho/s, theta, zeta) positions

Evaluate the basis for a list of (theta, zeta) positions on all flux surfaces given by s

Evaluate the basis and all derivatives for a list of (theta, zeta) positions on all flux surfaces given by s

Handle the selection of the functional and derivatives, based on the selection strings

Evaluate the basis with a tensorproduct for the given 1D (s, theta, zeta) values

Evaluate the basis with a tensorproduct for the given 1D (s, theta, zeta) values

Evaluate LA or NU and all derivatives for a list of (theta, zeta) positions on all flux surfaces given by s

Evaluate the mapping from reference to physical space (hmap)

Evaluate the mapping from reference to physical space (hmap) without logical coordinates

Evaluate the mapping from reference to physical space (hmap) without logical coordinates

Evaluate the mapping from reference to physical space (hmap)

evaluate the jacobian determinant and its derivatives

evaluate the jacobian determinant and its derivatives

evaluate components of the metric tensor and their derivatives

evaluate 1D-radial profiles and their derivatives with respect to rho

evaluate iota/pressure profile and its derivatives with respect to rho2=rho^2

Find which grid cell contains the given point

Find which grid cell contains the given point

Check if 1D spline belongs to span of given B-splines

Evaluate value of 1D spline at location x: y=S(x)

Evaluate derivative of 1D spline at location x: y=S'(x)

allocate all variables in fBase

change data from oldBase to self. Forier modes are directly copied so, if new mode space is smaller, its like a Fourier cut-off. if new modes do not match old ones, they are set to zero. Note that a change of nfp is not possible· as well as a change from sine to cosine

compare sf with the input type fBase

copy the type fBase

evaluate all modes at specific given point in theta and zeta

evaluate special 1D base in theta direction (cos(mt_i),sin(mt_i)) or its derivative(s) on a given set of points for tensor-product evaluation of 2D sin and cos base: sin(mthet-nzeta) = sin(mthet)cos(nzeta)-cos(mthet)sin(nzeta) == dot_product( (sin(mthet),-cos(mthet)) , (cos(nzeta),sin(nzeta))) cos(mthet-nzeta) = cos(mthet)cos(nzeta)+sin(mthet)sin(nzeta) == dot_product( (cos(mthet), sin(mthet)) , (cos(nzeta),sin(nzeta))) so for the 1D base, mTotal1d depends on using sin/cos/sin+cos base.

evaluate special 1D base in zeta direction (cos(mt_i),sin(mt_i)) or its derivative(s) on a given set of points for tensor-product evaluation of 2D sin and cos base: sin(mthet-nzeta) = sin(mthet)cos(nzeta)-cos(mthet)sin(nzeta) == dot_product( (sin(mthet),-cos(mthet)) , (cos(nzeta),sin(nzeta))) cos(mthet-nzeta) = cos(mthet)cos(nzeta)+sin(mthet)sin(nzeta) == dot_product( (cos(mthet), sin(mthet)) , (cos(nzeta),sin(nzeta))) so for the 1D base, nTotal1d is always 2*n_max+1

evaluate all modes at a list of given points in theta and zeta

evaluate all modes at all interpolation points

evaluate all modes at all interpolation points, making use of the tensor product: y_ij = DOFs_mn * SIN(mt_i - nz_j ) => SIN(mt_i) DOFs_mn COS(nz_j) -COS(mt_i) DOFs_mn SIN(nz_j) y_ij = DOFs_mn * COS(mt_i - nz_j ) => COS(mt_i) DOFs_mn COS(nz_j) +SIN(mt_i) DOFs_mn SIN(nz_j) => a1_im DOFs_mn b1_nj + a2_im DOFs_mn b2_nj can be written as 2 SPECIAL MATMAT operations: c(i,1,n)=a1(i,m) DOFs(m,n) , c(i,2,n) = a2(i,m) DOFs(m,n) => c(i,d,n) = DOT_PROD(a(i,d,1:mmax),DOFs(1:mmax,n)) y(i,j) = c(i,1,n) b1(n,j) + c(i,2,n) b2(n,j) = DOT_PROD(c(i,1:2,1:nmax),b(1:2,1:nmax,j) the 1D ordering in y does not neead a reshape, y(i,j) => y(1:mn_IP), 1D array data can be kept, as it is passed (with its start adress) to DGEMM.

evaluate all modes at a given interpolation point

evaluate all modes at a list of interpolation points

evaluate all modes on a tensor-produc grid (t_i,z_j), making use of the tensor product in the fourier series: y_ij = DOFs_mn * SIN(mt_i - nz_j ) => SIN(mt_i) DOFs_mn COS(nz_j) -COS(mt_i) DOFs_mn SIN(nz_j) y_ij = DOFs_mn * COS(mt_i - nz_j ) => COS(mt_i) DOFs_mn COS(nz_j) +SIN(mt_i) DOFs_mn SIN(nz_j) => a1_im DOFs_mn b1_nj + a2_im DOFs_mn b2_nj can be written as 2 SPECIAL MATMAT operations: c(i,1,n)=a1(i,m) DOFs(m,n) , c(i,2,n) = a2(i,m) DOFs(m,n) => c(i,d,n) = DOT_PROD(a(i,d,1:mmax),DOFs(1:mmax,n)) y(i,j) = c(i,1,n) b1(n,j) + c(i,2,n) b2(n,j) = DOT_PROD(c(i,1:2,1:nmax),b(1:2,1:nmax,j) the 1D ordering in y does not neead a reshape, y(i,j) => y(1:mn_IP), 1D array data can be kept, as it is passed (with its start adress) to DGEMM.

finalize the type fBase

initialize the type fBase maximum mode numbers, number of integration points, type of basis (sin/cos or sin and cos)

take values interpolated at sf%s_IP positions and project onto fourier basis by integration

allocate the type fBase

project from interpolation points to all modes DOFs = addDOFs+ fac MATMUL(base_IP_DOF,y_IP)

inverse of fBase_evalDOF_IP_tens

project from any 2D set of interpolation points, at tensor-product of (theta,zeta) positions given by "xn", to all modes DOFs = addDOFs+ fac MATMUL(base_xn,yn)

test fBase variable

Read ini file and put each line in a string object. All string objects are connected to a list of string objects starting with "firstString". MUST BE CALLED IN THE VERY BEGINNING OF THE PROGRAM!

Finalize Module

Finalize Module

Finalize Module

finalize the type functional

Finalize Module

Finalize Module

Finalize Module

Reset global variables

Finalize Module

Finalize VMEC module

on one flux surface, find for an given list of (thet_j,zeta_j), the corresponding (thet_j,zeta_j) positions, given Here, new boozer angles are theta=theta+Gt(theta,zeta) zeta=zeta+nu(theta,zeta), with Gt=lambda+iotanu and nu periodic functions and zero average and same base Note that in this routine, we will use a 2d root search with a newton method, setting [f1,f2]^T = [thet+A(thet,zeta)-thet=0, zeta+B(thet,zeta)-zeta*=0]^T that includes the derivatives (Jacobian), so that the newton step needs to the solved: -[f1] [ 1+dA/dthet dA/dzeta] [dthet] | | = | | | | -[f2] [ dB/dthet 1+dB/dzeta] [dzeta]

on one flux surface, find for a list of in thet_j,zeta_j, the corresponding (thet_j,zeta_j) positions, given Here, new PEST angles are theta=theta+lambda(theta,zeta) zeta=zeta, so a 1D root search in theta is is enough

Find parameter string containing keyword "Key" in list of strings starting with "FirstString" and return string "Str" without keyword. If keyword is not found in list of strings, return default values "Proposal" (error if not given). Ini file was read in before and is stored as list of character strings starting with "FirstString".

Fit disrete data along flux surfaces as spline for each fourier mode

Fit disrete data along flux surfaces as spline for each fourier mode input is given on the half mesh 2:nFluxVMEC

This routine was taken fom QUADRULE (http://people.sc.fsu.edu/~jburkardt/f_src/quadrule/quadrule.html)

This function returns the result of the 2D newton root search for the boozer angle

Derivative of the target function for finding the logical angle for given boozer angles

Target function for finding the logical angle for given boozer angles

Builds the boozer transform coordinate theta^B = theta + lambda + iota(s)*nu(s,theta,zeta) zeta^B = zeta +nu(s,theta,zeta)

Get command line arguments

compute different fields depending on the input parameters field_type and vector_component,

automatically build the string to be read from parameterfile, varname + m,n mode number, and then read it from parameterfile

Retrieve the integration points and weights (gauss points for radial integration)

This function returns the result of the 1D newton root search for the pest theta angle

Get unused file unit number

Read integer named "key" from setup file and store in "GETINT". If keyword "Key" is not found in ini file, the default value "Proposal" is used for "GETINT" (error if "Proposal" not given). Ini file was read in before and is stored as list of character strings starting with "FirstString".

Allocate and read integer array of unknown length "nIntegers" integer values named "Key" from ini file. If keyword "Key" is not found in setup file, the default values "Proposal" are used to create the array (error if "Proposal" not given). Setup file was read in before and is stored as list of character strings starting with "FirstString".

Read array of "nIntegers" integer values named "Key" from ini file. If keyword "Key" is not found in setup file, the default values "Proposal" are used to create the array (error if "Proposal" not given). Setup file was read in before and is stored as list of character strings starting with "FirstString".

Read logical named "key" from setup file and store in "GETINT". If keyword "Key" is not found in ini file, the default value "Proposal" is used for "GETINT" (error if "Proposal" not given). Ini file was read in before and is stored as list of character strings starting with "FirstString".

Return P L U matrices of the LU decomposition, cmoputed from LAPACK Routine (if P is not passed, L=P*L)

Create and initialize new string object.

Read real named "key" from setup file and store in "GETINT". If keyword "Key" is not found in ini file, the default value "Proposal" is used for "GETINT" (error if "Proposal" not given). Ini file was read in before and is stored as list of character strings starting with "FirstString".

Read array of "nReals" real values named "Key" from ini file. If keyword "Key" is not found in setup file, the default values "Proposal" are used to create the array (error if "Proposal" not given). Setup file was read in before and is stored as list of character strings starting with "FirstString".

Read array of "nReals" real values named "Key" from ini file. If keyword "Key" is not found in setup file, the default values "Proposal" are used to create the array (error if "Proposal" not given). Setup file was read in before and is stored as list of character strings starting with "FirstString".

Read string named "key" from setup file and store in "GETINT". If keyword "Key" is not found in ini file, the default value "Proposal" is used for "GETINT" (error if "Proposal" not given). Ini file was read in before and is stored as list of character strings starting with "FirstString".

Calculates current time (serial / OpenMP /MPI)

Calculates current time locally on a MPIrank (no MPI Barrier)

prepare all data to be written

write data to file

evaluate the mapping h (q1,q2,zeta) -> (x,y,z)

evaluate all metrics necessary for optimizer

evaluate all quantities at one given point (elemental)

evaluate the mapping h (q1,q2,zeta) -> (x,y,z)

evaluate total derivative of the mapping sum k=1,3 (dx(1:3)/dq^k) q_vec^k, where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) ,

evaluate total derivative of the mapping sum k=1,3 (dx(1:3)/dq^k) q_vec^k, where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) ,

evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,, where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,, where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate sum_k q_vec^k * sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) , k=1,2,3 where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate sum_k q_vec^k * sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) , k=1,2,3 where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(q^1,q^2,zeta)

evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(q^1,q^2,zeta)

evaluate derivative of Jacobian of mapping h: sum_k q_vec^k dJ_h/dq^k, k=1,2,3 at q=(q^1,q^2,zeta)

evaluate derivative of Jacobian of mapping h: sum_k q_vec^k dJ_h/dq^k, k=1,2,3 at q=(q^1,q^2,zeta)

evaluate curve X0(zeta), and T=X0',N,B,N',B', using the fourier series of X0_hat,N_hat and B_hat and transform from "hat" coordinates to cartesian coordinates: x=xhatcos(zeta)-sgn_rotyhatsin(zeta), y=yhatcos(zeta)+sgn_rotsin(zeta)xhat, z=zhat

finalize the type hmap_axisNB

evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 is evaluated at q_in=(X^1,X^2,zeta),

evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 is evaluated at q_in=(X^1,X^2,zeta),

evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta),

evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta),

initialize the type hmap_axisNB, reading from parameterfile and call init_params

initialize the aux variable

initialize the type hmap_axisNB and read "G-frame" from netcdf

test hmap_axisNB - evaluation of the map

evaluate the mapping h (X^1,X^2,zeta) -> (x,y,z)

evaluate all metrics necesseray for optimizer

evaluate all quantities at one given point (elemental)

evaluate total derivative of the mapping sum k=1,3 (dx(1:3)/dq^k) q_vec^k, where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) ,

evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,, where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) , k=1,2 where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(X^1,X^2,zeta)

evaluate derivative of Jacobian of mapping h: dJ_h/dq^k, k=1,2 at q=(X^1,X^2,zeta)

finalize the type hmap_cyl

evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 is evaluated at q_in=(X^1,X^2,zeta),

evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta),

initialize the type hmap_cyl, reading from parameterfile and then call init_params

Allocate and initialize auxiliary variable at zeta position.

initialize the type hmap_cyl, given the parameters as arguments

test hmap_cyl

evaluate the mapping h (X^1,X^2,zeta) -> (x,y,z) cartesian INFO: default routine that can be overwritten by specific hmap class, not using additional hmap-dependent auxiliary variables, but calling the pointwise routine eval

call %eval_aux on 1d array of points of size np, using auxiliary variable array of same size

evaluate total derivative of the mapping sum k=1,3 (dx(1:3)/dq^k) q_vec^k, where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) , INFO: default routine that can be overwritten by specific hmap class, not using additional hmap-dependent auxiliary variables, but calling the generic routine eval_dxdq_pw

call %eval_dxdq_aux on 1d array of points of size np, using auxiliary variable array of same size

evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,, where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta INFO: default routine that can be overwritten by specific hmap class, not using additional hmap-dependent auxiliary variables, but calling the pointwise routine eval_gij

evaluate sum_k sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) q_vec^k, k=1,2,3 where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta INFO: default routine that can be overwritten by specific hmap class, not using additional hmap-dependent auxiliary variables, but calling the pointwise routine eval_gij_dq

evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(X^1,X^2,zeta) INFO: default routine that can be overwritten by specific hmap class, not using additional hmap-dependent auxiliary variables, but calling the pointwise routine eval_Jh

call %eval_Jh_aux on 1d array of points of size np, using auxiliary variable array of same size

evaluate derivative of Jacobian of mapping h: sum_k dJ_h(q)/dq^k q_vec^k, k=1,2 at q=(X^1,X^2,zeta) INFO: default routine that can be overwritten by specific hmap class, not using additional hmap-dependent auxiliary variables, but calling the pointwise routine eval_Jh_dq

call %eval_Jh_dq1_aux on 1d array of points of size np, using auxiliary variable array of same size

evaluate the mapping h (q1,q2,zeta) -> (x,y,z)

evaluate all metrics necessary for optimizer

evaluate all quantities at one given point (elemental)

evaluate the mapping h (q1,q2,zeta) -> (x,y,z)

evaluate total derivative of the mapping sum k=1,3 (dx(1:3)/dq^k) q_vec^k, where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) ,

evaluate total derivative of the mapping sum k=1,3 (dx(1:3)/dq^k) q_vec^k, where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) ,

evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,, where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0 of ds and dtheta and 1 for dzeta

evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,, where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0 of ds and dtheta and 1 for dzeta

evaluate sum_k=1,3 q_vec^k * sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) , k=1,2 where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0 of ds and dtheta and 1 for dzeta

evaluate sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) , k=1,2 where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0 of ds and dtheta and 1 for dzeta

evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(q^1,q^2,zeta)

evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(q^1,q^2,zeta)

evaluate derivative of Jacobian of mapping h: sum_k q_vec^k * dJ_h/dq^k, k=1,2,3 at q=(q^1,q^2,zeta)

evaluate derivative of Jacobian of mapping h: sum_k q_vec^k * dJ_h/dq^k, k=1,2,3 at q=(q^1,q^2,zeta)

evaluate curve X0(zeta), position and first three derivatives, from given R0,Z0 Fourier

finalize the type hmap_frenet

evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 is evaluated at q_in=(X^1,X^2,zeta),

evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 is evaluated at q_in=(X^1,X^2,zeta),

evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta),

evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta),

initialize the type hmap_frenet, reading from parameterfile and then call init_params

initialize the aux variable

initialize the type hmap_frenet with number of elements

sign function depending on zeta, if omnig=False, sigma=1 if omnig=True, sigma=+1 for 0<=zeta<=pi/nfp, and -1 for pi/nfp<zeta<2pi

test hmap_frenet - evaluation of the map

evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 , at q_in=(X^1,X^2,zeta), INFO: default routine that can be overwritten by specific hmap class, not using additional hmap-dependent auxiliary variables, but calling the generic routine get_ddx_dqij

evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta), INFO: default routine that can be overwritten by specific hmap class, not using additional hmap-dependent auxiliary variables, but calling the generic routine get_dx_dqi

evaluate the mapping h (q1,q2,zeta) -> (x,y,z)

evaluate all metrics necesseray for optimizer

evaluate all quantities at one given point (elemental) NOTE: using calls to sf, not implemented/optimized for performance yet!

evaluate total derivative of the mapping sum k=1,3 (dx(1:3)/dq^k) q_vec^k, where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) ,

evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,, where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate sum_k sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) *q_vec^k, k=1,2,3 where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(q^1,q^2,zeta)

evaluate derivative of Jacobian of mapping h: sum_k q_vec^k * dJ_h/dq^k, k=1,2,3 at q=(q^1,q^2,zeta)

evaluate the effective major radius coordinate Rl(q)

evaluate the effective vertical coordinate Zl(q)

finalize the type hmap_knot

evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 is evaluated at q_in=(X^1,X^2,zeta),

evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta),

initialize the type hmap_knot, reading from parameter file and then call init_params

Allocate and initialize auxiliary variable at zeta position.

initialize the type hmap_knot, from given parameters as arguments

test hmap_knot - evaluation of the map

initialize the type hmap, also readin parameters here if necessary

initialize the hmap auxiliary variables, depends on hmap type

evaluate the mapping h (X^1,X^2,zeta) -> (x,y,z) cartesian

evaluate all metrics necesseray for optimizer

evaluate all quantities at one given point (elemental)

evaluate the mapping h (X^1,X^2,zeta) -> (x,y,z) cartesian INFO: overwrites default routine hmap_eval_aux from c_hmap.f90

evaluate total derivative of the mapping sum k=1,3 (dx(1:3)/dq^k) q_vec^k, where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) ,

evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,, where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate sum_k sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) *q_vec^k, k=1,2,3 where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and dzeta_dalpha then known to be either 0.0 for ds and dtheta and 1.0 for dzeta

evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(X^1,X^2,zeta)

evaluate derivative of Jacobian of mapping h: sum q_vec^k * dJ_h/dq^k, k=1,2,3 at q=(X^1,X^2,zeta)

finalize the type hmap_RZ

evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 is evaluated at q_in=(X^1,X^2,zeta),

evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta),

initialize the type hmap_RZ, no additional readin from parameterfile needed.

Allocate and initialize auxiliary variable at zeta position.

test hmap_RZ

Prints out remaining strings in list after read-in is complete

initialize base classes declared in _vars module, needed for computation of output fields

initialize a SFL-Boozer object, with some parameters taken from the state (globals) Note: as of v0.2.16 (Thanks Christopher Albert) f90wrap supports ALLOCATABLEs in the return value

Initialize Module

Initialize LAMBDA FROM U_init%X1,%X2 and iota profile, this computation is distributed over MPIranks

Initialize Module

Overwrite axis with average axis by center of closed line of the boundary in each poloidal plane

initialize the type functional with number of elements

Initialize Module

Build a 1D Vandermonde matrix using the Lagrange basis functions of degree N_In, evaluated at the interpolation points xi_Out

Initialize Module

Initialize Module

Initialise Profiles at GP!!!

Initialize Module

Initialize the solution with the given boundary condition

Initialize Module

Initialize VMEC module

used as rhs, then overwritten in solve

Computes matrix inverse using LAPACK Input matrix should be a square matrix

compute inverse of 3x3 matrix, needs sDet=1/det(Mat)

Computes all Lagrange functions evaluated at position x in [-1,1] For details see paper Barycentric Lagrange Interpolation by Berrut and Trefethen (SIAM 2004) Uses function ALMOSTEQUAL Algorithm 34, Kopriva book

Solve for lambda on one given flux surface (spos_in), using weak form of J^s=0: d/dzeta(B_theta)-d/dtheta(B_zeta)=0

Compute Legendre-Gauss nodes and integration weights (algorithm 23, Kopriva book)

Evaluate the Legendre polynomial L_N and its derivative at position x[-1,1] recursive algorithm using the N_in-1 N_in-2 Legendre polynomials algorithm 22, Kopriva book

Starting with initial guess by Parter Relation, a Newton method is used to find the roots of the Legendre Polynomial Lder_(N_in), which are the positions of Gauss-Lobatto points. Uses qAndLEvaluation subroutine. algorithm 25, Kopriva

Transform upper case letters in "Str1" into lower case letters, result is "Str2", but only up the the equal sign.

Compute Equilibrium, iteratively

Compute Equilibrium, iteratively

Computes mth polynomial differentiation matrix for interpolation polynomial given by set of nodes. (Algorithm 38, Kopriva book)

close netcdf file

define a dimension to the netCDF file

after creating a new file and making all definitions, one has to call end_def_mode

if variable name contains "/", these are interpreted as groups/subgroups. split the varname at first occurence of "/" to get the first group name on the file level. Then get the group id. repeat until no "/" is found anymore. output the final groupid and the variable name without the group names.

closes file and frees variable

get integer or real array of dimension 1d,2d,3d,4d (depends on optional argument) netcdf call get_var knows type and dimensions directly from argument abort if variable does not exist. USE var_exists for checking

get integer or real scalar (depends on optional argument) abort if variable does not exist. USE var_exists for checking

get the size of a ulti-dimensional array for all dimensions ndims

get the number of dimensions of a variable

netcdf error handling via sf%ioError variable

allocate and initialize class and open/close the netcdf file and define read ("r") or write ("w" includes read) mode

open netcdf file

define and put an array value to the netCDF file

define and put a char to the netCDF file as a string NOTE: This is only used for naming coordinate directions with a single character

define and put a scalar value to the netCDF file

check if variable name exists (also including groups separated with "/")

d/dx f(x) function of the Min1D type wrapped for Root1D

f(x) function of the Min1D type wrapped for Root1D

constructor for the Min1D type wrapped for Root1D

Newton's iterative algorithm to find the minimimum of function f(x) in the interval [a,b], using df(x)=0 and the derivative

Newton's iterative algorithm to find the minimimum of function f(x,y) in the interval x(i)[a(i),b(i)], using grad(f(x)=0 and the derivative

Newton's iterative algorithm to find the root of function FR(x(:)) in the interval [a(:),b(:)], using d/dx(:)F(x)=0 and the derivative

Newton's iterative algorithm to find the root of function FR(x(:)) in the interval [a(:),b(:)], using d/dx(:)F(x)=0 and the derivative

Newton's iterative algorithm to find the root of function [f1(x1,x2),f2(x1,x2)]=[0,0] in the interval a(i)<=x(i)<=b(i), using the Jacobian dfi/dxj, i=1,2, j=1,2, such that fi(x1,x2)=fi(x1_0,x2_0)+ [dfi/dx1,dfi/dx2].[dx1,dx2] in each step, we find dx1,dx2 st -[[dfi/dxj]] dxj =fi(x1_0,x2_0)

normalizes a nDim vector with respect to the eucledian norm

Find MAX/MIN/SUM of 1D array (assumed-shape) across all MPI ranks and bradcast result back to all MPI ranks.

Find MAX/MIN/SUM of 2D array (assumed-shape) across all MPI ranks and bradcast result back to all MPI ranks.

Find MAX/MIN/SUM scalar value across MPI ranks and bradcast result back to all MPI ranks.

Find MAX/MIN/SUM scalar value across MPI ranks and bradcast result back to all MPI ranks.

Barrier for specified communicator, or world-communicator otherwise.

Broadcast a 1D array (assumed-shape) from MPI rank 'fromRank' to all MPI ranks.

Broadcast a 1D array (assumed-shape) from MPI rank 'fromRank' to all MPI ranks.

Broadcast a 1D array (assumed-shape) from MPI rank 'fromRank' to all MPI ranks.

Broadcast a 2D array (assumed-shape) from MPI rank 'fromRank' to all MPI ranks.

Broadcast a scalar from MPI rank 'fromRank' to all MPI ranks.

Broadcast a scalar from MPI rank 'fromRank' to all MPI ranks.

Broadcast a scalar from MPI rank 'fromRank' to all MPI ranks.

Deinitialisation of MPI.

Broadcast a 1D array (assumed-shape) from MPI rank 'fromRank' to all MPI ranks (nonblocking)

Broadcast a 2D array (assumed-shape) from MPI rank 'fromRank' to all MPI ranks (nonblocking)

Initialisation of MPI.

Find on MPI rank 'toRank' MAX/MIN/SUM of 1D array (assumed-shape) across all MPI ranks (nonblocking).

Find on MPI rank 'toRank' MAX/MIN/SUM of 2D array (assumed-shape) across all MPI ranks.

Find on MPI rank 'toRank' MAX/MIN/SUM of 1D array (assumed-shape) across all MPI ranks.

Find on MPI rank 'toRank' MAX/MIN/SUM of 2D array (assumed-shape) across all MPI ranks.

Find on MPI rank 'toRank' MAX/MIN/SUM scalar value across MPI ranks.

Find on MPI rank 'toRank' MAX/MIN/SUM scalar value across MPI ranks.

Wait for completion of a single nonblocking communication

Wait for completion of all nonblocking communications for req(:)

Sum an array across MPI ranks: explicit-shape with implicit reshaping Multi-D->1D.

Function for 1D newton root search for PEST

calculate the prefactor for the d-th coefficient of the n-th derivative of a polynomial

Computes polynomial differentiation matrix for interpolation polynomial given by set of nodes. (Algorithm 37, Kopriva book)

get the exact polynomial antiderivative, with respect to rho2

evaluate the n-th derivative of a power polynomial

finalize the type rProfile

Print a progress bar to screen, call either with init=T or init=F

Evaluate the polynomial q=L_{N_in+1}-L_{N_in-1} and its derivative at position x in [-1,1] Recursive algorithm using the N_in-1 N_in-2 Legendre polynomials. (Algorithm 24, Kopriva book)

READ axis from netcdf file, needs netcdf library! ======= HEADER OF THE NETCDF FILE VERSION 3.1 =================================================================================== === FILE DESCRIPTION: * axis, normal and binormal of the frame are given in cartesian coordinates along the curve parameter zeta [0,2pi]. * The curve is allowed to have a field periodicity NFP, but the curve must be provided on a full turn. * The adata is given in real space, sampled along equidistant zeta point positions: zeta(i)=(i+0.5)/nzeta * (2pi/NFP), i=0,...,nzeta-1 always shifted by (2pi/NFP) for the next field period. Thus the number of points along the axis for a full turn is NFP*nzeta * definition of the axis-following frame in cartesian coordinates ( boundary surface at rho=1):

READ FROM SPECIFIC NETCDF FILE: general data and "boundary" group ======= HEADER OF THE NETCDF FILE VERSION 3.1 =================================================================================== === FILE DESCRIPTION: * axis, normal and binormal of the frame are given in cartesian coordinates along the curve parameter zeta [0,2pi]. * The curve is allowed to have a field periodicity NFP, but the curve must be provided on a full turn. * The data is given in REAL SPACE, sampled along equidistant zeta point positions: zeta(i)=(i+0.5)/nzeta * (2pi/NFP), i=0,...,nzeta-1 always shifted by (2pi/NFP) for the next field period. Thus the number of points along the axis for a full turn is NFP*nzeta * definition of the axis-following frame in cartesian coordinates ( boundary surface at rho=1):

read an input solution and initialize U(0) (X1,X2,LA) of size X1/X2/LA_base , from an ascii .dat file if size of grid/X1/X2/LA not equal X1/X2/X3_base interpolate readin solution to the current base of Uin

read an input solution and initialize U(0) (X1,X2,LA) of size X1/X2/LA_base , from an ascii .dat file if size of grid/X1/X2/LA not equal X1/X2/X3_base interpolate readin solution to the current base of Uin

evaluate the n-th derivative of (rho^2) with respect to rho ~sqrt(magnetic flux).

Rodrigues' rotation formula assumption for now is that the origin is fixed at rot_origin=(/0.,0.,0./) and the rotation axis is fixed at rot_axis=(/0.,0.,1./)

evaluate the 2nd derivative of a radial profile with respect to rho ~sqrt(magnetic flux).

evaluate the 3rd derivative of a radial profile with respect to rho ~sqrt(magnetic flux).

evaluate the 4th derivative of a radial profile with respect to rho ~sqrt(magnetic flux)

evaluate the n-th derivative of a radial profile with respect to rho ~sqrt(magnetic flux). NOTE: n has to be in [0,4] due to an explicit implementation of the product rule.

Evaluate value at x of all basis functions with support in local cell values[j] = B_j(x) for jmin <= j <= jmin+degree

Evaluate value and n derivatives at x of all basis functions with support in local cell derivs[i,j] = (d/dx)^i B_j(x) for 0 <= i <= n and jmin <= j <= jmin+degree

Evaluate derivative at x of all basis functions with support in local cell derivs[j] = B_j'(x) for jmin <= j <= jmin+degree

Free storage

Initialize non-uniform B-splines object

Evaluate value at x of all basis functions with support in local cell values[j] = B_j(x) for jmin <= j <= jmin+degree

Evaluate value and n derivatives at x of all basis functions with support in local cell derivs[i,j] = (d/dx)^i B_j(x) for 0 <= i <= n and jmin <= j <= jmin+degree

Evaluate derivative at x of all basis functions with support in local cell derivs[j] = B_j'(x) for jmin <= j <= jmin+degree

Free storage

Initialize uniform B-splines object

Private subroutine for assembling and factorizing linear system for any kind of boundary conditions at xmin and xmax

Evaluate value of 1D spline at all locations in array x

Evaluate derivative of 1D spline at all locations in array x

Destroy 1D spline (re-initialization is possible afterwards)

Initialize 1D spline object as element of span(B-splines)

Compute interpolating 1D spline values on grid. If Hermite BCs are used, function derivatives at appropriate boundary are also needed.

Destroy local objects and free allocated memory

Get coordinates of interpolation points (1D grid)

Initialize a 1D spline interpolator object

allocate all variables in sbase

apply boundary conditions at axis and edge, via solving the Lagrange multiplier problem: x_new=x_old & A*x_new = d

apply strong boundary conditions at axis and edge Not used anymore, WAS FOUND TO BE NOT STABLE (OSCILLATORY) IN GENERAL, especially for BCs with only derivatives prescribed... new LGM method below is used now!

apply strong boundary conditions at axis and edge for solution update

change data from old_sBase to self. using interpolations of the old data at the new interpolation points

compare sf with input sbase

copy onto sf <-- tocopy

evaluate sbase at position x [0,1], NOT EFFICIENT!!

simply evaluate function or derivative at point x, for multiple DOF vectors

simply evaluate function with a base or base derivative evaluated at a point and its corresponding iElem use together with sBase_eval(x) => iElem,base

evaluate all degrees of freedom at all Gauss Points (deriv=0 solution, deriv=1 first derivative d/ds)

simply evaluate function or derivative at point x

finalize the type sBase

initialize the type sbase with polynomial degree, continuity ( -1: disc, 1: full) and number of gauss points per element

take values interpolated at sf%s_IP positions and give back the degrees of freedom

initialize the type sbase with polynomial degree, continuity ( -1: disc, 1: full) and number of gauss points per element

test sbase variable

Handle the selection of the base, based on the selection string

Handle the selection of the base, based on the selection string

interface to find_boozer_angles from the class t_sfl_boozer

interface to find_boozer_angles from the class t_sfl_boozer

finalize sfl boozer class

initialize sfl boozer class

compare to sf grid with input grid to see if they are the same

copy the type sgrid, copies sf <= tocopy ... call sf%copy(tocopy)

find grid cell for certain position

finalize the type sgrid

initialize the type sgrid with number of elements

test sgrid variable

Allocate and initialize uniform or non-uniform B-splines

@brief Print error message to standard-error, stop execution and dump backtrace information. @details

Calculate number of cells from number of interpolation points calculate the number of grid cells that yields the desired number of interpolation points

@brief Print warning message to standard-error and continue execution. @details

res=aX+bY , where X,Y,res are of type t_var_sol

copy tocopy => sf

free (=deallocate) sf of type t_sol_var_MHD3D

initialize (=allocate) sf of type t_sol_var

|X|^2, where X is of type t_var_sol, so three values are returned: |X1|^2,|X2|^2,|LA|^2

set all variables to scalar

set variabes X1,X2,LA of toset => sf, optional argument to scale toset with a scalar (for example -1.0_wp)

test sol_var_MHD3D

Solve linear system of dimension dims and multiple RHS

Solve linear system of dimension dims and multiple RHS

Helper routine to go from spline interpolation points to spline coefficients and apply smooth axis boundary condition.

Transform a function from the GVEC angles q(s,theta,zeta) to new angles q(s,theta,zeta*) by using interpolation in angular direction (fourier transform) and spline interpolation in radial direction (at s_IP points of output base) the interpolation points are given by thetazeta_IP, which are the angle positions of an equidistant interpolation grid in PEST/Boozer angles

Transform a function from VMEC angles q(s,theta,zeta) to new angles q(s,theta,zeta) by projection onto the modes of the new angles: sigma_mn(theta,zeta) using a given in s Here, new angles are theta=theta+A(theta,zeta), zeta=zeta+B(theta,zeta), with A,B periodic functions and zero average and same base Note that in this routine, the integral is transformed back to (theta,zeta) q_mn = iint_0^2pi q(theta,zeta) sigma_mn(theta,zeta) dtheta dzeta = iint_0^2pi q(theta,zeta) sigma_mn(theta,zeta) [(1+dA/dtheta)(1+dB/dzeta)-(dA/dzetadB/dzeta)] dtheta dzeta

get_new

Allocate class and call init

Write evaluation of the axis and signed axisNB frame to file

Write evaluation of the axis and signed frenet frame to file

Visualize VMEC flux surface data for each mode, for Rmnc

Visualize VMEC flux surface data in planes or 3D, number of radial posisiotns fixed to nFluxVMEC+1, only R,Z,Lambda

evaluate spline for specific mode at position s

write modes prepared above

Subroutine to write

Subroutine to write multidimensional data to netCDF format

Subroutine to write 3D point data to VTK format

write an input solution (X1,X2,LA) to an ascii .dat file