SPEC 3.20
Stepped Pressure Equilibrium Code
Functions/Subroutines
Diagnostics to check the code

More...

Functions/Subroutines

subroutine bfield (zeta, st, Bst)
 Compute the magnetic field. More...
 
subroutine hesian (NGdof, position, Mvol, mn, LGdof)
 Computes eigenvalues and eigenvectors of derivative matrix, \(\nabla_{\bf xi}{\bf F}\). More...
 
subroutine jo00aa (lvol, Ntz, lquad, mn)
 Measures error in Beltrami equation, \(\nabla \times {\bf B} - \mu {\bf B}\). More...
 
subroutine pp00aa
 Constructs Poincaré plot and "approximate" rotational-transform (driver). More...
 
subroutine pp00ab (lvol, sti, Nz, nPpts, poincaredata, fittedtransform, utflag)
 Constructs Poincaré plot and "approximate" rotational-transform (for single field line). More...
 
subroutine stzxyz (lvol, stz, RpZ)
 Calculates coordinates, \({\bf x}(s,\theta,\zeta) \equiv R \, {\bf e}_R + Z \, {\bf e}_Z\), and metrics, at given \((s,\theta,\zeta)\). More...
 

Detailed Description

Function/Subroutine Documentation

◆ bfield()

subroutine bfield ( real, intent(in)  zeta,
real, dimension(1:node), intent(in)  st,
real, dimension(1:node), intent(out)  Bst 
)

Compute the magnetic field.

Returns the magnetic field field line equations, \( d{\bf x}/d\phi = {\bf B}/B^\phi \) .

Equations of field line flow

  • The equations for the fieldlines are normalized to the toroidal field, i.e.

    \begin{equation} \dot s \equiv \frac{B^s }{B^\zeta}, \qquad \dot \theta \equiv \frac{B^\theta}{B^\zeta}. \label{eq:stdot_bfield} \end{equation}

Representation of magnetic field

  • The components of the vector potential, \({\bf A}=A_\theta \nabla + A_\zeta \nabla \zeta\), are

    \begin{eqnarray} A_\theta(s,\theta,\zeta) &=& \sum_{i,l} {\color{red} A_{\theta,e,i,l}} \; {\overline T}_{l,i}(s) \cos\alpha_i + \sum_{i,l} {\color{Orange} A_{\theta,o,i,l}} \; {\overline T}_{l,i}(s) \sin\alpha_i, \label{eq:At_bfield} \\ A_\zeta( s,\theta,\zeta) &=& \sum_{i,l} {\color{blue} A_{\zeta, e,i,l}} \; {\overline T}_{l,i}(s) \cos\alpha_i + \sum_{i,l} {\color{Cerulean}A_{\zeta ,o,i,l}} \; {\overline T}_{l,i}(s) \sin\alpha_i, \label{eq:Az_bfield} \end{eqnarray}

    where \({\overline T}_{l,i}(s) \equiv \bar s^{m_i/2} \, T_l(s)\), \(T_l(s)\) is the Chebyshev polynomial, and \(\alpha_j \equiv m_j\theta-n_j\zeta\). The regularity factor, \(\bar s^{m_i/2}\), where \(\bar s \equiv (1+s)/2\), is only included if there is a coordinate singularity in the domain (i.e. only in the innermost volume, and only in cylindrical and toroidal geometry.)
  • The magnetic field, \(\sqrt g \, {\bf B} = \sqrt g B^s {\bf e}_s + \sqrt g B^\theta {\bf e}_\theta + \sqrt g B^\zeta {\bf e}_\zeta\), is

    \begin{eqnarray} \begin{array}{ccccrcrcrcrcccccccccccccccccccccccccccccccccccccccccccccccccccc} \sqrt g \, {\bf B} & = & {\bf e}_s & \sum_{i,l} [ ( & - m_i {\color{blue} A_{\zeta, e,i,l}} & - & n_i {\color{red} A_{\theta,e,i,l}} & ) {\overline T}_{l,i} \sin\alpha_i + ( & + m_i {\color{Cerulean}A_{\zeta ,o,i,l}} & + & n_i {\color{Orange} A_{\theta,o,i,l}} & ) {\overline T}_{l,i} \cos\alpha_i ] \\ & + & {\bf e}_\theta & \sum_{i,l} [ ( & & - & {\color{blue} A_{\zeta, e,i,l}} & ) {\overline T}_{l,i}^\prime \cos\alpha_i + ( & & - & {\color{Cerulean}A_{\zeta ,o,i,l}} & ) {\overline T}_{l,i}^\prime \sin\alpha_i ] \\ & + & {\bf e}_\zeta & \sum_{i,l} [ ( & {\color{red} A_{\theta,e,i,l}} & & & ) {\overline T}_{l,i}^\prime \cos\alpha_i + ( & {\color{Orange} A_{\theta,o,i,l}} & & & ) {\overline T}_{l,i}^\prime \sin\alpha_i ] \end{array} \end{eqnarray}

  • In Eqn. \((\ref{eq:stdot_bfield})\) , the coordinate Jacobian, \(\sqrt g\), cancels. No coordinate metric information is required to construct the fieldline equations from the magnetic vector potential.

IT IS REQUIRED TO SET IVOL THROUGH GLOBAL MEMORY BEFORE CALLING BFIELD.

The format of this subroutine is constrained by the NAG ode integration routines.

Parameters
[in]zetatoroidal angle \( \zeta \)
[in]stradial coordinate \(s\) and poloidal angle \(\theta\)
[out]Bsttangential magnetic field directions \(B_s, B_\theta\)

References allglobal::ate, allglobal::ato, allglobal::aze, allglobal::azo, bfield(), allglobal::cpus, allglobal::gbzeta, get_cheby(), get_zernike(), constants::half, allglobal::halfmm, allglobal::im, allglobal::in, allglobal::ivol, allglobal::lcoordinatesingularity, inputlist::lrad, allglobal::mn, allglobal::mpi_comm_spec, inputlist::mpol, allglobal::myid, allglobal::ncpu, allglobal::node, allglobal::notstellsym, constants::one, fileunits::ounit, allglobal::regumm, numerical::small, constants::two, numerical::vsmall, inputlist::wmacros, and constants::zero.

Referenced by bfield(), bfield_tangent(), jo00aa(), pp00ab(), and sphdf5::write_grid().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ hesian()

subroutine hesian ( integer, intent(in)  NGdof,
real, dimension(0:ngdof)  position,
integer, intent(in)  Mvol,
integer, intent(in)  mn,
integer, intent(in)  LGdof 
)

Computes eigenvalues and eigenvectors of derivative matrix, \(\nabla_{\bf xi}{\bf F}\).

Parameters
[in]NGdofnumber of global degrees of freedom
[in,out]positioninternal geometrical degrees of freedom
[in]Mvoltotal number of volumes in computation
[in]mnnumber of Fourier harmonics
[in]LGdofwhat is this?

construction of Hessian matrix

  • The routine dforce() is used to compute the derivatives, with respect to interface geometry, of the force imbalance harmonics, \([[p+B^2/2]]_{j}\), which may be considered to be the "physical" constraints, and if Igeometry==3 then also the derivatives of the "artificial" spectral constraints, \(I_j \equiv (R_\theta X + Z_\theta Y)_j\).
  • The input variable Lconstraint determines how the enclosed fluxes, \(\Delta \psi_t\) and \(\Delta \psi_p\), and the helicity multiplier, \(\mu\), vary as the geometry is varied; see global.f90 and mp00ac() for more details.

construction of eigenvalues and eigenvectors

  • If LHevalues==T then the eigenvalues of the Hessian are computed using the NAG routine F02EBF.
  • If LHevectors==T then the eigenvalues and the eigenvectors of the Hessian are computed.
  • Note that if Igeometry==3, then the derivative-matrix also contains information regarding how the "artificial" spectral constraints vary with geometry; so, the eigenvalues and eigenvectors are not purely "physical".

  • The eigenvalues and eigenvectors (if required) are written to the file .ext.GF.ev as follows:

    open(hunit,file="."//trim(ext)//".GF.ev",status="unknown",form="unformatted")
    write(hunit)ngdof,ldvr,ldvi ! integers; if only the eigenvalues were computed then Ldvr=Ldvi=1;
    write(hunit)evalr(1:ngdof) ! reals ; real part of eigenvalues;
    write(hunit)evali(1:ngdof) ! reals ; imaginary part of eigenvalues;
    write(hunit)evecr(1:ngdof,1:ngdof) ! reals ; real part of eigenvalues; only if Ldvr=NGdof;
    write(hunit)eveci(1:ngdof,1:ngdof) ! reals ; imaginary part of eigenvalues; only if Ldvi=NGdof;
    close(hunit)
  • The eigenvectors are saved in columns of evecr, eveci, as described by the NAG documentation for F02EBF.

References allglobal::cpus, allglobal::dbbdmp, allglobal::dbbdrz, allglobal::denergydrr, allglobal::denergydrz, allglobal::denergydzr, allglobal::denergydzz, allglobal::dessian, allglobal::dessian2d, allglobal::dffdrz, dforce(), allglobal::dmupfdx, inputlist::dpp, inputlist::dqq, allglobal::drbc, allglobal::drbs, allglobal::dzbc, allglobal::dzbs, allglobal::energy, constants::half, allglobal::hdffdrz, inputlist::helicity, hesian(), allglobal::hessian, allglobal::hessian2d, fileunits::hunit, inputlist::igeometry, allglobal::im, allglobal::in, allglobal::irbc, allglobal::irbs, allglobal::izbc, allglobal::izbs, allglobal::lbbintegral, inputlist::lcheck, inputlist::lfindzero, inputlist::lfreebound, allglobal::lhessian2dallocated, allglobal::lhessian3dallocated, allglobal::lhessianallocated, inputlist::lhevalues, inputlist::lhevectors, inputlist::lhmatrix, allglobal::localconstraint, inputlist::lperturbed, allglobal::mpi_comm_spec, inputlist::mu, fileunits::munit, allglobal::myid, allglobal::ncpu, allglobal::notstellsym, inputlist::nvol, constants::one, fileunits::ounit, packxi(), inputlist::pflux, preset(), allglobal::psifactor, numerical::small, numerical::sqrtmachprec, constants::ten, constants::two, numerical::vsmall, inputlist::wmacros, allglobal::yesstellsym, and constants::zero.

Referenced by hesian(), and spec().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ jo00aa()

subroutine jo00aa ( integer, intent(in)  lvol,
integer, intent(in)  Ntz,
integer, intent(in)  lquad,
integer, intent(in)  mn 
)

Measures error in Beltrami equation, \(\nabla \times {\bf B} - \mu {\bf B}\).

This routine is called by xspech() as a post diagnostic and only if Lcheck==1.

construction of current, \({\bf j} \equiv \nabla \times \nabla \times {\bf A}\)

  • The components of the vector potential, \({\bf A}=A_\theta \nabla + A_\zeta \nabla \zeta\), are

    \begin{eqnarray} A_\theta(s,\theta,\zeta) &=& \sum_{i,l} {\color{red} A_{\theta,e,i,l}} \; {\overline T}_{l,i}(s) \cos\alpha_i + \sum_{i,l} {\color{Orange} A_{\theta,o,i,l}} \; {\overline T}_{l,i}(s) \sin\alpha_i, \label{eq:At_jo00aa} \\ A_\zeta (s,\theta,\zeta) &=& \sum_{i,l} {\color{blue} A_{\zeta, e,i,l}} \; {\overline T}_{l,i}(s) \cos\alpha_i + \sum_{i,l} {\color{Cerulean}A_{\zeta ,o,i,l}} \; {\overline T}_{l,i}(s) \sin\alpha_i, \label{eq:Az_jo00aa} \end{eqnarray}

    where \({\overline T}_{l,i}(s) \equiv \bar s^{m_i/2} \, T_l(s)\), \(T_l(s)\) is the Chebyshev polynomial, and \(\alpha_j \equiv m_j\theta-n_j\zeta\). The regularity factor, \(\bar s^{m_i/2}\), where \(\bar s \equiv (1+s)/2\), is only included if there is a coordinate singularity in the domain (i.e. only in the innermost volume, and only in cylindrical and toroidal geometry.)
  • The magnetic field, \(\sqrt g \, {\bf B} = \sqrt g B^s {\bf e}_s + \sqrt g B^\theta {\bf e}_\theta + \sqrt g B^\zeta {\bf e}_\zeta\), is

    \begin{eqnarray} \begin{array}{ccccrcrcrcrcccccccccccccccccccccccccccccccccccccccccccccccccccc} \sqrt g \, {\bf B} & = & {\bf e}_s & \sum_{i,l} [ ( & - m_i {\color{blue} A_{\zeta, e,i,l}} & - & n_i {\color{red} A_{\theta,e,i,l}} & ) {\overline T}_{l,i} \sin\alpha_i + ( & + m_i {\color{Cerulean}A_{\zeta ,o,i,l}} & + & n_i {\color{Orange} A_{\theta,o,i,l}} & ) {\overline T}_{l,i} \cos\alpha_i ] \\ & + & {\bf e}_\theta & \sum_{i,l} [ ( & & - & {\color{blue} A_{\zeta, e,i,l}} & ) {\overline T}_{l,i}^\prime \cos\alpha_i + ( & & - & {\color{Cerulean}A_{\zeta ,o,i,l}} & ) {\overline T}_{l,i}^\prime \sin\alpha_i ] \\ & + & {\bf e}_\zeta & \sum_{i,l} [ ( & {\color{red} A_{\theta,e,i,l}} & & & ) {\overline T}_{l,i}^\prime \cos\alpha_i + ( & {\color{Orange} A_{\theta,o,i,l}} & & & ) {\overline T}_{l,i}^\prime \sin\alpha_i ] \end{array} \end{eqnarray}

  • The current is

    \begin{eqnarray} \sqrt g \, {\bf j} = ( \partial_\theta B_\zeta - \partial_\zeta B_\theta) \; {\bf e}_s + ( \partial_\zeta B_s - \partial_s B_\zeta ) \; {\bf e}_\theta + ( \partial_s B_\theta - \partial_\theta B_s ) \; {\bf e}_\zeta , \end{eqnarray}

    where (for computational convenience) the covariant components of \({\bf B}\) are computed as

    \begin{eqnarray} B_s & = & (\sqrt g B^s) \, g_{ss } / \sqrt g + (\sqrt g B^\theta) \, g_{s \theta} / \sqrt g + (\sqrt g B^\zeta) \, g_{s \zeta} / \sqrt g, \\ B_\theta & = & (\sqrt g B^s) \, g_{s\theta} / \sqrt g + (\sqrt g B^\theta) \, g_{\theta\theta} / \sqrt g + (\sqrt g B^\zeta) \, g_{\theta\zeta} / \sqrt g, \\ B_\zeta & = & (\sqrt g B^s) \, g_{s\zeta } / \sqrt g + (\sqrt g B^\theta) \, g_{\theta\zeta } / \sqrt g + (\sqrt g B^\zeta) \, g_{\zeta \zeta} / \sqrt g. \end{eqnarray}

quantification of the error

  • The measures of the error are

    \begin{eqnarray} ||\left( {\bf j}-\mu {\bf B}\right)\cdot\nabla s || & \equiv & \int \!\! ds \oint\!\!\!\oint \!d\theta d\zeta \,\,\, \left| \sqrt g \, {\bf j}\cdot\nabla s -\mu \; \sqrt g \, {\bf B} \cdot\nabla s \right|, \label{eq:Es_jo00aa} \\ ||\left( {\bf j}-\mu {\bf B}\right)\cdot\nabla \theta || & \equiv & \int \!\! ds \oint\!\!\!\oint \!d\theta d\zeta \,\,\, \left| \sqrt g \, {\bf j}\cdot\nabla \theta-\mu \; \sqrt g \, {\bf B} \cdot\nabla \theta \right|, \label{eq:Et_jo00aa} \\ ||\left( {\bf j}-\mu {\bf B}\right)\cdot\nabla \zeta || & \equiv & \int \!\! ds \oint\!\!\!\oint \!d\theta d\zeta \,\,\, \left| \sqrt g \, {\bf j}\cdot\nabla \zeta -\mu \; \sqrt g \, {\bf B} \cdot\nabla \zeta \right|. \label{eq:Ez_jo00aa} \end{eqnarray}

comments

  • Is there a better definition and quantification of the error? For example, should we employ an error measure that is dimensionless?
  • If the coordinate singularity is in the domain, then \(|\nabla\theta|\rightarrow\infty\) at the coordinate origin. What then happens to \(||\left( {\bf j}-\mu {\bf B}\right)\cdot\nabla \theta ||\) as defined in Eqn. \((\ref{eq:Et_jo00aa})\)?
  • What is the predicted scaling of the error in the Chebyshev-Fourier representation scale with numerical resolution? Note that the predicted error scaling for \(E^s\), \(E^\theta\) and \(E^\zeta\) may not be standard, as various radial derivatives are taken to compute the components of \({\bf j}\). (See for example the discussion in Sec.IV.C in Hudson et al. (2011) [4] , where the expected scaling of the error for a finite-element implementation is confirmed numerically.)
  • Instead of using Gaussian integration to compute the integral over \(s\), an adaptive quadrature algorithm may be preferable.
Parameters
[in]lvolin which volume should the Beltrami error be computed
[in]Ntznumber of grid points in \(\theta\) and \(\zeta\)
[in]lquaddegree of Gaussian quadrature
[in]mnnumber of Fourier harmonics

details of the numerics

  • The integration over \(s\) is performed using Gaussian integration, e.g., \(\displaystyle \int \!\! f(s) ds \approx \sum_k \omega_k f(s_k)\); with the abscissae, \(s_k\), and the weights, \(\omega_k\), for \(k=1,\) Iquad \(_v\), determined by CDGQF. The resolution, N \( \equiv \) Iquad \(_v\), is determined by Nquad (see global.f90 and preset() ). A fatal error is enforced by jo00aa() if CDGQF returns an ifail \(\ne 0\).

  • Inside the Gaussian quadrature loop, i.e. for each \(s_k\),

    • The metric elements, \(g_{\mu,\nu} \equiv \) gij(1:6,0,1:Ntz), and the Jacobian, \(\sqrt g \equiv \) sg(0,1:Ntz), are calculated on a regular angular grid, \((\theta_i,\zeta_j)\), in coords(). The derivatives \(\partial_i g_{\mu,\nu} \equiv\) gij(1:6,i,1:Ntz) and \(\partial_i \sqrt g \equiv \) sg(i,1:Ntz), with respect to \( i \in \{ s,\theta,\zeta \}\) are also returned.

    • The Fourier components of the vector potential given in Eqn. \((\ref{eq:At_jo00aa})\) and Eqn. \((\ref{eq:Az_jo00aa})\), and their first and second radial derivatives, are summed.

    • The quantities \(\sqrt g B^s\), \(\sqrt g B^\theta\) and \(\sqrt g B^\zeta\), and their first and second derivatives with respect to \((s,\theta,\zeta)\), are computed on the regular angular grid (using FFTs).

    • The following quantities are then computed on the regular angular grid

      \begin{eqnarray} \sqrt g j^s & = & \sum_u \left[ \partial_\theta (\sqrt g B^u) \; g_{u,\zeta } + (\sqrt g B^u) \; \partial_\theta g_{u,\zeta } - (\sqrt g B^u) g_{u,\zeta } \; \partial_\theta \sqrt g / \sqrt g \right] / \sqrt g \nonumber \\ & - & \sum_u \left[ \partial_\zeta (\sqrt g B^u) \; g_{u,\theta} + (\sqrt g B^u) \; \partial_\zeta g_{u,\theta} - (\sqrt g B^u) g_{u,\theta} \; \partial_\zeta \sqrt g / \sqrt g \right] / \sqrt g, \\ \sqrt g j^\theta & = & \sum_u \left[ \partial_\zeta (\sqrt g B^u) \; g_{u,s } + (\sqrt g B^u) \; \partial_\zeta g_{u,s } - (\sqrt g B^u) g_{u,s } \; \partial_\zeta \sqrt g / \sqrt g \right] / \sqrt g \nonumber \\ & - & \sum_u \left[ \partial_s (\sqrt g B^u) \; g_{u,\zeta } + (\sqrt g B^u) \; \partial_s g_{u,\zeta } - (\sqrt g B^u) g_{u,\zeta } \; \partial_s \sqrt g / \sqrt g \right] / \sqrt g, \\ \sqrt g j^\zeta & = & \sum_u \left[ \partial_s (\sqrt g B^u) \; g_{u,\theta} + (\sqrt g B^u) \; \partial_s g_{u,\theta} - (\sqrt g B^u) g_{u,\theta} \; \partial_s \sqrt g / \sqrt g \right] / \sqrt g \nonumber \\ & - & \sum_u \left[ \partial_\theta (\sqrt g B^u) \; g_{u,s } + (\sqrt g B^u) \; \partial_\theta g_{u,s } - (\sqrt g B^u) g_{u,s } \; \partial_\theta \sqrt g / \sqrt g \right] / \sqrt g. \end{eqnarray}

  • The error is stored into an array called beltramierror which is then written to the HDF5 file in hdfint().

References allglobal::ate, allglobal::ato, allglobal::aze, allglobal::azo, allglobal::beltramierror, bfield(), allglobal::cfmn, allglobal::cheby, coords(), allglobal::cpus, allglobal::dpflux, allglobal::dtflux, allglobal::efmn, allglobal::gbzeta, get_cheby_d2(), get_zernike_d2(), allglobal::guvij, constants::half, inputlist::igeometry, allglobal::im, allglobal::in, allglobal::ivol, jo00aa(), allglobal::lcoordinatesingularity, inputlist::lerrortype, inputlist::lrad, allglobal::mpi_comm_spec, inputlist::mpol, inputlist::mu, allglobal::myid, inputlist::nfp, allglobal::node, allglobal::notstellsym, allglobal::nt, inputlist::nvol, allglobal::nz, allglobal::ofmn, constants::one, fileunits::ounit, constants::pi2, allglobal::regumm, allglobal::rij, allglobal::rtt, allglobal::sfmn, allglobal::sg, allglobal::tt, constants::two, inputlist::wmacros, allglobal::zernike, constants::zero, and allglobal::zij.

Referenced by final_diagnostics(), and jo00aa().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ pp00aa()

subroutine pp00aa

Constructs Poincaré plot and "approximate" rotational-transform (driver).

relevant input variables

  • The resolution of Poincaré plot is controlled by
    • nPtraj trajectories will be located in each volume;
    • nPpts iterations per trajectory;
    • odetol o.d.e. integration tolerance;
  • The magnetic field is given by bfield() .
  • The approximate rotational transform is determined, in pp00ab() , by fieldline integration.

format of output: Poincaré

  • The Poincaré data is written to .ext.poincare:xxxx , where xxxx is an integer indicating the volume. The format of this file is as follows:

    write(svol,'(i4.4)')lvol ! lvol labels volume;
    open(lunit+myid,file="."//trim(ext)//".poincare."//svol,status="unknown",form="unformatted")
    do until end of file
    write(lunit+myid) nz, nppts ! integers
    write(lunit+myid) data(1:4,0:nz-1,1:nppts) ! doubles
    enddo
    close(lunit+myid)

    where

    • \(\theta \equiv\,\)data(1,k,j) is the poloidal angle,
    • \( s \equiv\,\)data(2,k,j) is the radial coordinate,
    • \( R \equiv\,\)data(3,k,j) is the cylindrical \(R\),
    • \( Z \equiv\,\)data(4,k,j) is the cylindrical \(Z\),
  • The integer k=0,Nz-1 labels toroidal planes, so that \(\phi = ( 2 \pi / \texttt{Nfp}) ( k / \texttt{Nz})\),
  • The integer j=1,nPpts labels toroidal iterations.
  • Usually (if no fieldline integration errors are encountered) the number of fieldlines followed in volume lvol is given by \(N+1\), where the radial resolution, \(N\equiv\,\)Ni(lvol) , is given on input. This will be over-ruled by if nPtrj(lvol) , given on input, is non-negative.
  • The starting location for the fieldline integrations are equally spaced in the radial coordinate \(s_i=s_{l-1}+ i (s_{l}-s_{l-1})/N\) for \(i=0,N\), along the line \(\theta=0\), \(\zeta=0\).

format of output: rotational-transform

  • The rotational-transform data is written to .ext.transform:xxxx , where xxxx is an integer indicating the volume. The format of this file is as follows:
    open(lunit+myid,file="."//trim(ext)//".sp.t."//svol,status="unknown",form="unformatted")
    write(lunit+myid) lnptrj-ioff+1 ! integer
    write(lunit+myid) diotadxup(0:1,0,lvol) ! doubles
    write(lunit+myid) ( fiota(itrj,1:2), itrj = ioff, lnptrj ) ! doubles
    close(lunit+myid)

References allglobal::cpus, allglobal::diotadxup, sphdf5::finalize_flt_output(), constants::half, inputlist::igeometry, sphdf5::init_flt_output(), inputlist::iota, allglobal::ivol, inputlist::lconstraint, allglobal::lcoordinatesingularity, allglobal::lmns, allglobal::lplasmaregion, inputlist::lrad, allglobal::lvacuumregion, allglobal::mpi_comm_spec, allglobal::myid, allglobal::ncpu, inputlist::nppts, inputlist::nptrj, inputlist::nvol, allglobal::nz, inputlist::odetol, inputlist::oita, constants::one, fileunits::ounit, constants::pi, pp00aa(), pp00ab(), inputlist::ppts, constants::two, inputlist::wmacros, sphdf5::write_poincare(), sphdf5::write_transform(), and constants::zero.

Referenced by final_diagnostics(), and pp00aa().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ pp00ab()

subroutine pp00ab ( integer, intent(in)  lvol,
real, dimension(1:2)  sti,
integer, intent(in)  Nz,
integer, intent(in)  nPpts,
real, dimension(1:4,0:nz-1,1:nppts)  poincaredata,
real, dimension(1:2)  fittedtransform,
integer, intent(out)  utflag 
)

Constructs Poincaré plot and "approximate" rotational-transform (for single field line).

relevant input variables

  • The resolution of Poincaré plot is controlled by
    • nPpts iterations per trajectory;
    • odetol o.d.e. integration tolerance;
    The magnetic field is given by bfield() .

rotational-transform

  • The approximate rotational transform is determined by field line integration. This is constructed by fitting a least squares fit to the field line trajectory.
Parameters
[in]lvol
sti
[in]Nz
[in]nPpts
poincaredata
fittedtransform
[out]utflag

References bfield(), allglobal::cpus, allglobal::ivol, allglobal::mpi_comm_spec, allglobal::myid, allglobal::ncpu, allglobal::node, inputlist::nvol, inputlist::odetol, constants::one, fileunits::ounit, constants::pi2, pp00ab(), numerical::small, stzxyz(), constants::two, volume(), and constants::zero.

Referenced by pp00aa(), and pp00ab().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ stzxyz()

subroutine stzxyz ( integer, intent(in)  lvol,
real, dimension(1:3), intent(in)  stz,
real, dimension(1:3), intent(out)  RpZ 
)

Calculates coordinates, \({\bf x}(s,\theta,\zeta) \equiv R \, {\bf e}_R + Z \, {\bf e}_Z\), and metrics, at given \((s,\theta,\zeta)\).

  • This routine is a "copy" of co01aa(), which calculates the coordinate information on a regular, discrete grid in \(\theta\) and \(\zeta\) at given \(s\) whereas stzxyz() calculates the coordinate information at a single point \((s,\theta,\zeta)\).
  • Todo:
    Please see co01aa() for documentation.

Parameters
[in]lvol
[in]stz
[out]RpZ

References allglobal::cpus, constants::half, allglobal::halfmm, inputlist::igeometry, allglobal::im, allglobal::in, allglobal::irbc, allglobal::irbs, allglobal::izbc, allglobal::izbs, allglobal::lcoordinatesingularity, allglobal::mn, allglobal::mpi_comm_spec, allglobal::myid, allglobal::notstellsym, inputlist::ntor, inputlist::nvol, constants::one, fileunits::ounit, stzxyz(), numerical::vsmall, and constants::zero.

Referenced by pp00ab(), and stzxyz().

Here is the call graph for this function:
Here is the caller graph for this function: