Geometry

Functions/Subroutines
subroutine	coords (lvol, lss, Lcurvature, Ntz, mn)
	Calculates coordinates, \({\bf x}(s,\theta,\zeta) \equiv R \, {\bf e}_R + Z \, {\bf e}_Z\), and metrics, using FFTs. More...

Detailed Description

Function/Subroutine Documentation

coords()

subroutine coords	(	integer, intent(in)	lvol,
		real, intent(in)	lss,
		integer, intent(in), value	Lcurvature,
		integer, intent(in)	Ntz,
		integer, intent(in)	mn
	)

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

Coordinates

• We work in coordinates, \((s,\theta,\zeta)\), which are be defined inversely via a transformation to Cartesian coordinates, \((x,y,z)\).
• The toroidal angle, \(\zeta\), is identical to the cylindrical angle, \(\zeta\equiv\phi\).
• The radial coordinate, \(s\), is not a global variable: it only needs to be defined in each volume, and in each volume \(s \in [-1,1]\).
• The choice of poloidal angle, \(\theta\), does not affect the following.

Geometry

• The geometry of the "ideal"-interfaces, \({\bf x}_v(\theta,\zeta)\), is given by \(R(\theta,\zeta)\) and \(Z(\theta,\zeta)\) as follows:

  • Igeometry=1 : Cartesian

    \begin{eqnarray} {\bf x} & \equiv & r_{pol}\theta \; {\bf \hat i} + r_{tor}\zeta \; {\bf \hat j}+ R \; {\bf \hat k} \end{eqnarray}

    where \(r_{pol}\) and \(r_{tor}\) are inputs and \(r_{pol}=r_{tor}=1\) by default.
  • Igeometry=2 : Cylindrical

    \begin{eqnarray} {\bf x} & = & R \; \cos\theta \; {\bf \hat i} + R \; \sin\theta \; {\bf \hat j} + \zeta \; {\bf \hat k} \end{eqnarray}

  • Igeometry=3 : Toroidal

    \begin{eqnarray} {\bf x} & \equiv & R \; {\bf \hat r} + Z \; {\bf \hat k} \end{eqnarray}

    where \({\bf \hat r} \equiv \cos \phi \; {\bf \hat i} + \sin \phi \; {\bf \hat j}\) and \({\bf \hat \phi} \equiv - \sin \phi \; {\bf \hat i} + \cos \phi \; {\bf \hat j}\).

• The geometry of the ideal interfaces is given as Fourier summation: e.g., for stellarator-symmetry

  \begin{eqnarray} R_v(\theta,\zeta) & \equiv & \sum_j R_{j,v} \cos\alpha_j, \\ Z_v(\theta,\zeta) & \equiv & \sum_j Z_{j,v} \sin\alpha_j, \end{eqnarray}

  where \(\alpha_j \equiv m_j \theta - n_j \zeta\).

interpolation between interfaces

• The "coordinate" functions, \(R(s,\theta,\zeta)\) and \(Z(s,\theta,\zeta)\), are constructed by radially interpolating the Fourier representations of the ideal-interfaces.
• The \(v\)-th volume is bounded by \({\bf x}_{v-1}\) and \({\bf x}_{v}\).
• In each annular volume, the coordinates are constructed by linear interpolation:

  \begin{eqnarray} \begin{array}{cccccccccccccccccccccccccccc} R(s,\theta,\zeta) & \equiv & \displaystyle \sum_j & \displaystyle \left[ \; \frac{(1-s)}{2} \; R_{j,v-1} + \frac{(1+s)}{2} \; R_{j,v}\; \right] \; \cos\alpha_j,\\ Z(s,\theta,\zeta) & \equiv & \displaystyle \sum_j & \displaystyle \left[ \; \frac{(1-s)}{2} \; Z_{j,v-1} + \frac{(1+s)}{2} \; Z_{j,v}\; \right] \; \sin\alpha_j, \end{array} \end{eqnarray}

coordinate singularity: regularized extrapolation

• For cylindrical or toroidal geometry, in the innermost, "simple-torus" volume, the coordinates are constructed by an interpolation that "encourages" the interpolated coordinate surfaces to not intersect.
• Introduce \(\bar s \equiv (s+1)/2\), so that in each volume \(\bar s \in [0,1]\), then

  \begin{eqnarray} R_{j}(s) & = & R_{j,0} + (R_{j,1} - R_{j,0} ) f_j, \\ Z_{j}(s) & = & Z_{j,0} + (Z_{j,1} - Z_{j,0} ) f_j, \end{eqnarray}

  where, in toroidal geometry,

  \begin{eqnarray} f_j \equiv \left\{ \begin{array}{llcccccccccccccc} \bar s & , & \mbox{ for } m_j=0, \\ \bar s^{m_j} & , & \mbox{ otherwise.} \end{array}\right\}. \end{eqnarray}

• Note: The location of the coordinate axis, i.e. the \(R_{j,0}\) and \(Z_{j,0}\), is set in the coordinate "packing" and "unpacking" routine, packxi().

Jacobian

• The coordinate Jacobian (and some other metric information) is given by

  • Igeometry=1 : Cartesian

    \begin{eqnarray} {\bf e}_\theta \times {\bf e}_\zeta & = & -r_{tor}R_\theta \; \hat {\bf i} - r_{pol}R_\zeta \; \hat {\bf j} + r_{pol}r_{tor}\hat {\bf k} \\ \boldsymbol{\xi} \cdot {\bf e}_\theta \times {\bf e}_\zeta & = & \delta R \\ \sqrt g & = & R_s \ r_{pol} \ r_{tor} \end{eqnarray}

  • Igeometry=2 : Cylindrical

    \begin{eqnarray} {\bf e}_\theta \times {\bf e}_\zeta & = & (R_\theta \sin \theta + R \cos\theta ) \; {\bf \hat i} + (R \sin \theta - R_\theta \cos\theta ) \; {\bf \hat j} - R R_\zeta \; {\bf \hat k} \\ \boldsymbol{\xi}\cdot {\bf e}_\theta \times {\bf e}_\zeta & = & \delta R \; R \\ \sqrt g & = & R_s \; R \end{eqnarray}

  • Igeometry=3 : Toroidal

    \begin{eqnarray} {\bf e}_\theta \times {\bf e}_\zeta & = & - R \, Z_\theta \, \hat r + (Z_\theta \,R_\zeta - R_\theta \,Z_\zeta) \hat \phi + R \,R_\theta \,\hat z\\ \boldsymbol{\xi}\cdot {\bf e}_\theta \times {\bf e}_\zeta & = & R ( \delta Z \; R_\theta - \delta R \; Z_\theta ) \\ \sqrt g & = & R ( Z_s \; R_\theta - R_s \; Z_\theta ) \end{eqnarray}

cartesian metrics

• The cartesian metrics are

  \begin{eqnarray} \begin{array}{cccccccccccccccccccccccccccccccccccccccccccccccc} g_{s s } = R_s R_s ,& g_{s \theta} = R_s R_\theta ,& g_{s \zeta } = R_s R_\zeta ,& g_{\theta\theta} = R_\theta R_\theta + r_{pol}^2 ,& g_{\theta\zeta } = R_\theta R_\zeta ,& g_{\zeta \zeta } = R_\zeta R_\zeta + r_{tor}^2 \end{array} \end{eqnarray}

cylindrical metrics

• The cylindrical metrics are

  \begin{eqnarray} \begin{array}{cccccccccccccccccccccccccccccccccccccccccccccccc} g_{s s } = R_s R_s ,& g_{s \theta} = R_s R_\theta ,& g_{s \zeta } = R_s R_\zeta ,& g_{\theta\theta} = R_\theta R_\theta + R^2 ,& g_{\theta\zeta } = R_\theta R_\zeta ,& g_{\zeta \zeta } = R_\zeta R_\zeta + 1 \end{array} \end{eqnarray}

logical control

• The logical control is provided by Lcurvature as follows:
  • Lcurvature=0 : only the coordinate transformation is computed, i.e. only \(R\) and \(Z\) are calculated, e.g. global()
  • Lcurvature=1 : the Jacobian, \(\sqrt g \), and "lower" metrics, \(g_{\mu,\nu}\), are calculated , e.g. bnorml(), lforce(), curent(), metrix(), sc00aa()
  • Lcurvature=2 : the "curvature" terms are calculated, by which I mean the second derivatives of the position vector; this information is required for computing the current, \({\bf j}=\nabla\times\nabla\times{\bf A}\), e.g. jo00aa()
  • Lcurvature=3 : the derivative of the \(g_{\mu,\nu}/\sqrt g\) w.r.t. the interface boundary geometry is calculated, e.g. metrix(), curent()
  • Lcurvature=4 : the derivative of the \(g_{\mu,\nu}\) w.r.t. the interface boundary geometry is calculated, e.g. dforce()
  • Lcurvature=5 : the derivative of \(\sqrt{g}\) w.r.t. the interface boundary geometry is calculated, e.g. rzaxis()

Parameters

[in]	lvol	specified in which volume to compute coordinates
[in]	lss	radial coordinate \(s\)
[in]	Lcurvature	logical control flag
[in]	Ntz	number of points in \(\theta\) and \(\zeta\)
[in]	mn	number of Fourier harmonics

References coords(), allglobal::cosi, allglobal::cpus, allglobal::dbdx, allglobal::drodr, allglobal::drodz, allglobal::dzodr, allglobal::dzodz, allglobal::guvij, constants::half, allglobal::halfmm, inputlist::igeometry, allglobal::im, allglobal::in, allglobal::irbc, allglobal::irbs, allglobal::izbc, allglobal::izbs, allglobal::lcoordinatesingularity, allglobal::mpi_comm_spec, allglobal::myid, allglobal::notstellsym, allglobal::nt, inputlist::ntor, allglobal::nz, constants::one, fileunits::ounit, constants::pi2, allglobal::rij, inputlist::rpol, inputlist::rtor, allglobal::sg, allglobal::sini, numerical::small, constants::two, volume(), numerical::vsmall, inputlist::zbc, inputlist::zbs, constants::zero, and allglobal::zij.

Referenced by compute_guvijsave(), coords(), curent(), jo00aa(), lforce(), preset(), rzaxis(), volume(), and sphdf5::write_grid().

Here is the call graph for this function:

Here is the caller graph for this function: